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ABSTRACT 


This report documents the Zero-g Thermodynamic Venting System (TVS) performance 
prediction computer program. The zero-g TVS is a device that destratifies and rejects 
environmentally induced zero-g thermal gradients in the LH 2 storage transfer system. A 
recirculation pump and spray injection manifold recirculates liquid throughout the length of the 
tank, thereby destratifying both the ullage gas and liquid bulk. Heat rejection is accomplished by 
the opening of the TVS control valve which allows a small flow rate to expand to a low pressure 
thereby producing a low temperature heat sink which is used to absorb heat from the recirculating 
liquid flow. The program was written in FORTRAN 77 language on the HP-9000 and IBM PC 
computers. It can be run on various platforms with a FORTRAN compiler. 
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SECTION 1 
INTRODUCTION 


1.1 Purpose 

The purpose of the zero-g thermodynamic venting system (TVS) model is to define the pessure 
level requirements, propellant loss due to venting, total pump power consumption, venting system 
operation duration and frequency as a function of liquid level and acceleration environment. 


1.2 Problem Description 

Long-term storage of subcritical cryogens in space is subjected to thermal stratification which is 
more severe than experienced in a 1-g environment due to the absence of gravity-induced body 
forces. If left uncontrolled, the thermal gradients result in excessive tank pressure rise and 
formation of liquid/vapor mixtures within the liquid bulk, Hquid acquisition device, and propellant 
transfer lines. A subsystem is th^fore needed to reduce the therm^ gradient to acceptable levels, 
and reject the environmental heat leakage in an efficient manner. 


1.3 Areas of Application 

This program is designed to predict the tank ullage pressure and tenq)erature, pump flow rate and 
I^wer consumption, propellant loss due to venting, venting duration and frequrency for different 
liquid levels, accelerations and heat leak rates. It can also be used to noodel tank chill down and 
liquid fill in zero-g and one-g environments. The progrpi therefore has general applicability in 
areas of long-term cryogenic fluid storage and transfer in space and on the ground. 


1.4 Description of the Physical Model 

The zero-g TVS concept, shown in Figure 1.1, has been defined to operate in a self-induced, 
forced convection environment It includes a recirculation pump, located external to the tank, to 
flow liquid fiom the tank, a spray manifold and injection tube to mix and destratify the ullage and 
liquid bulk, a heat exchanger to abscnrb heat fiom the tank, and an overboard vent line to reject the 
heat. The concept operates in the following manner. When the vent pressure level is reached, the 
recirculation pump is activated, resulting in liquid flow through the spray injection manifold which 
destratifies the liquid bulk and ullage gas through forced convection. V^en the fluid bulk 
temperature reaches a predetermine level, the vent valve is opened, resulting in a temperature drop 
of the vent flow by isenthalpic expansion. This liquid is used to absorb heat fiom the recirculation 
flow via the heat exchanger, and subsequently vented overboard as a gas. 
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SECTION 2 

ANALYTICAL MODEL DESCRIPTION 


The zero-g TVS model consists of the thermal-fluid models of the heat exchanger, spray manifold 
and injection tubes, recirculation pump, and tank. These models were developed and verified 
independently before they were integrated into the transient TVS model. Following is a description 
of each model. 


2.1 Heat Exchanger Model 

The heat exchanger model is based on the generalized two-phase cryogenic propellant dump model 
developed to evaluate the Space Shuttle Main Propulsion System (MPS) cryogenic propellant 
dump/vacuum inerting operations performance (Ref. 29). It is a multi-node finite difference model 
that simulates two-phase flow in a quasi steady-state mode. 

The model uses the fluid properties at the inlet of the spray manifold as the input to the first node. 
With one of the inlet fluid property, namely the total enthalpy of the fluid at the inlet, the total 
enthalpy at the exit can be calculate based on the First Law of Thermodynamics. The total 
enthalpy at the exit, and an assumed mass flow rate are used to determine the exit static pressure. 
The exit static pressure is determined by an iterative process: with the flow assum^ choked at the 
exit, the exit static pressure is increased incrementally until the maximum entropy is achieved 
(sonic flow), or when it becomes grwter than the back pressure (subsonic flow). From the 
calculated exit pressure, the other exit fluid properties of the last node, the total pressure loss ^ 
between the inlet and outlet can then be calculated and the inlet fluid properties can be determined. 

The following sections will provide the equations used in the heat exchanger model. 

2.1.1 Huid Quality At Heat Exchanger Outlet 

The outlet static pressure is calculated, assuming choked or sonic flow, where the entropy point is 
mayimiim. The following equations arc solved simultaneously for the liquid quality of the fluid at 
the outlet 


hQ = hj + ~ + AHa 
m 


Vo = 


m 

PoA 


Po = 


(Pl)< 


+ Y, 


(Pv)< 


1 

1 

(Pl)o. 


*>o = (1 - yoX*>L)o + Yo(hv)o 

where the following fluid properties ate based on the outlet static pressure 


( 2 . 1 . 1 ) 

( 2 . 1 . 2 ) 

(2.1.3) 


(2.1.4) 
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(hL)Q= the outlet liquid enthalpy 
(h V )q = the outlet vapor enthalpy 
(Pl)o= the outlet liquid density 
the outlet vapor density 

Q = total heat-transfer rate to a specific node 

AH= change in height of the line between inlet and outlet 

g 

a = — = acceleration 

gc 

ho= total enthalpy at the outlet 
Po- total density at the oudet 
Vo=fluid velocity at the outlet 
Yq= fluid quality at the outlet 

With the oudet quality, the total entropy then can be calculated using the following equation 
{(S)o = (l-yoXSL)o + Yo(Sv)o}_ (2.1.5) 


where (Sl)q= the oudet liquid entropy 
(Sv)o= tile oudet vapor entropy 

It^tion of the above oudet equations can be performed to obtain the maximum entropy point and 
the oudet static pressure. 

2.1.2 Two-Phase Pressure Loss in the Heat Exchang er 

To calculate the two-phase pressure loss (momentum and fiiction) between the iidet and oudet of 
the heat exchanger, the Lockhart-Martinelli correlation is used. The oudet pressure is 

Po=(Ps)o + (Pd)o (2.1.6) 

The total pressure loss term is further defined as 

APt=AP„ + AP, (2.1.7) 


where AP„ ^pressure loss due to momentum change 
APf = pressure loss due to frictional forces 

The momentum pressure loss is defined as 

AP„=^(Vo-V,) (2.1.8) 


The fiictional pressure loss is defined as 


Ift'' ' 


m 


IP: 


ft 

I 
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144K 

1 

1 

2Pl8c 

A 


<I>l' 


where K = “ j = the line loss coefficient 

Pl= the average liquid density between the inlet and outlet 

Y = the average total liquid quality between the inlet and outlet 
d>L=f(X) is the Lockhart-Martinelli correlation factor 

The Lockhart-Martinelli correlation is approximately defined as 

^ X 
X is defined as 

0.5 




0.1 


0.9 

f- \ 

X = 



1-Y 


Pv 




1 Y J 


^Pl> 


where Py = the average vapor densiQ^ between the inlet and outlet 
the average liquid viscosity between the inlet and outlet 
p.y = the average vapor viscosity between the inlet and outlet 

2.1.3 Forced Convection Heat -Transfer Model 

The heat-transfer equations used in the steady-state model are as follows 

Two-phase heat transfer using the correlation proposed by John C. Chen (1963) 

J = [hp<.F + h^S]AT 


where 


hpc= 0.0231 




k, j UJ 


^ 0.»p 0.45„ 0.49 OJS . j0.24^0.« 
hp2 - 0.00122—*' ^0.3p^^0.29j^0.24p^0.2* 


(2.1.9) 


(2.1.10) 

( 2 . 1 . 11 ) 


( 2 . 1 . 12 ) 

(2.1.13) 

(2.1.14) 
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F = f(X„) 


(2.1.15) 



(2.1.16) 

M-l 

AT = T*-Ts 

(2.1.17) 

A^p _ ATpyX, 
T. 

(2.1.18) 


s 


The single-phase heat-transfer coirelation used in the mcxlel (liquid and superheated gas) 



(2.1.19) 



-2.1.4 Spray Temperature 


( 2 . 1 . 20 ) 


To shorten the run time of the zero-g TVS model, an alternative was provided to calculate the spray 
temperature as a function of the liquid subcooling and tank pressure. TTiis is obtained from an 
energy balance between the hot and cold fluids of the heat exchanger 


Ts=Tp- 


mv Ah 


mp c„. 


( 2 . 1 . 21 ) 


where Tj is the spray temperature 
Tp is the pump temperature 

mv is the vent flow rate 
mp is the pump flow rate 

Ah is the heat absorption capability of the vent flow 
Cp^ is the heat ctpacity of the liquid 

The TVS vent flow rate and heat absorption c^ability were calculated as a function of the liquid 
subcooling and tank flow rate and are shown in Figs. 2.1 and 2.2, respectively. Tliis data is 
provided as a table look-up to the zero-g TVS model. 

2.2 Spray Manifold and Injection Tube Model 

Fluid is recirculated from the tank to the spoy manifold and injection tubes where it is sprayed into 
the ullage and liquid. A one-dimensional, incompressible fluid dynamic model was developed to 
determine the pressures in the spray manifold and injection tubes, and to calculate the spray flow 
rates and velocities leaving the injection orifices. Following is a description of the model 
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LIQUID SUBCOOLING - (psid) 


— 1 
1 0 


( 

Hgurc 2. 1 .2 Heat Absorption Gipability of Vent Flow as a Function of Liquid Subcooling 
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2.2.1 


The spray manifold calculates the pressure drop through the manifold and determines the pressure 
at the inlet of the spray injection mbes (Fig. 2.2.1). The model accounts for frictional pressure 
drop in the manifold, and pressure losses resulting from flow turning and contraction at the exit of 
the inanifold. From Bernoulli equation 

+ + + ( 2 . 2 . 1 ) 

P P ^8c 

where (PsM)i the spray manifold inlet pressure 

(Psm)o the spray manifold outlet pressure 
Vs,^, is the velocity in the spray manifold 
Zi, Zq are the inlet and outlet elevations 

a = — is the acceleration 

gc 

The total head loss is defined as 


(^l)sm”^sm 2g 


(2.2.2) 

The total loss coefficient Ks,^( is given by 


Ksm (^f )sM 

SM 

(2.2.3) 

and includes 




(spray manifold frictional loss coefficient ) 

(2.2.4a) 


(90-degree bend resistance at the manifold exit ) 

(2.2.4b) 

C=.)„-a=[i-(| 

^ (sudden contraction at the manifold exit ) 

mJ ^ 

(2.2.4c) 

In Eqs. 2.2.4, 

^SM “ 

is the spray manifold length 



is the bend equivalent length 
Djj is the spray injection tube ID 
Dsm is the spray manifold ID 

f is the friction coefficient in the spray manifold obtained frx>m 


n-7 




1 


Re >3000 



(2.2.5) 


f.„ = — ,ReS3000 


Re 


Eq. 2.2. 1 can be solved for the spray manifold outlet pressure 
(Psm)o “ (PsM)i “ ^SM^SM ” P^SM 

where the dynamic pressure in the spray manifold is given by 


V 2 1 

Qsm = P^ = ^ 

2gc 2pg, 




m« 


^Asm^ 


(2.2.6) 


(2.2.7) 


2.2.2 Spray Injection Tube 

The spray injection tube model is a multinode noodel which assigns a node to each orifice (Fig 
2.2.2). Bernoulli equation is first applied to find the pressure do^stream of the inlet 90 degree 
bend of the injection tube (pressure at die inlet of the straight section) 


Pi (Psm)o ^i(^b)si 


( 2 . 2 . 8 ) 


In Eq. 2.2.8, is the 90 degree bend resistance and q^ is the inlet dynamic pressure given by 


qi = 


1 

2pgc 


r • V 


JEl 



(2.2.9) 


where Aj, is the flow area of an injection tube and m^ is the mass flow rate in each tube (equal to 
the flow rate in the manifold divided by the number of tubes). 


The straight secticHi of the spray injection tube is divided into 45 equal inodes corresponding to the 
45 spray o^ces. Each node has a pressure and a mass flow rate at the inlet (i), cen^, and outlet 
(o) of the node. The outlet pressure and mass flow rate of one node is therefore the inlet pressure 
and mass flow rate of the preceding node 


Bernoulli equation is applied successively from inlet to center, and from center to outlet to 
determine the pressure at the center and outlet of a node il 
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From inlet to center, 

Pn=(Pi)„+paY-K,(q,)„ 

where Az is the nodal length and Kf is the frictional loss coefficient. 
From center to outlet, 

(Po), = P.+paY-Kf(q,)„ 


where the outlet dynamic pressure of node n is given by 



1 

(“»1 

2pg. 

Asi 


( 2 . 2 . 10 ) 


( 2 . 2 . 11 ) 


( 2 . 2 . 12 ) 


n-10 




The mass flow rate at node n outlet [ m^ ] is obtained from 



(2.2.13) 


^mg j in Eq. 2.2.13 is the spray flow rate calculated from an incompressible flow relation 


(ms) =(As), 


1 


(2pg.[p.-(PT)J 


K. 


(2.2.14) 


In Eq. 2.2.14, Kg is the loss coefficient of an orifice in a duct given by 



(2.2.15) 


where is the discharge coefficient ( Q ==0.8), and As/At is the ratio of the orifice to the tank 
area (As/At= 0). Thus, K, is determined to be 1.56. 


The tank pressure (p^ at node n is calculated as 
(Pt )a = Pu (ullage nodes) 

= Pu + pLg^a, (liquid nodes) (2.2. 16) 

where is the distance from the liquid surface to node n. 

2.2.3 Spray Manif old and Injection Tube Model Algorithm 

The flow chart of the spray manifold and injection tube model is given in Section 3.1.2. The 
model starts out with a guess of the pump flow rate and calculates the pressures and mass flow 
rates at each node. Knowing the pressure and spray flow rate of the last node N, it then calculates 
the tank pressure corresponding to that last node by solving the inccMiipressible flow relation of Eq. 



(2.2.17) 


Next, (Pt)n.«ic compared with (pt)^^ obtained from the ullage pressure and hydrostatic head 
(Eq. 2.2.16). If they are not equal witl^ a specified tolerance (0.001 psi), a new guess of the 
pump flow rate will be made and the process repeated until convergence on (pt)j^ is achieved. 


2.3 Recirculation Pump Model 
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The zero-g TVS LH2 recirculation pump is a centrifugal pump which is a constant output pressure 
device since it imparts kinetic pressure to the fluid due to rotation. Cons^uently, the pump 
pressure rise ( Ap^) is only a function of rotation speed (N) and tip velocity (U) 

7^ (2.3.1) 

720 

where is the impeller diameter. 

The fluid horsepower required by the pump flow ( m ), raised to Ap^ pressure, is equal to 


^ mApp (2.3.2) 

° flpP 

where qp is the pump mechanical efficiency. 

The pump operating speed then changes as a result of the energy absorbed by the fluid and the 
power supplied to the pump through a power source. The instantaneous rate of change m pump 
operating speed is 


dN f HPt-HPp ^ 
dt I, IpN , 


6.0185x10^ 


(2.3.3) 


where Ip is the polar nooment of inertia of the pump and HPj is die input power to the pump. 
Integration of the pump acceleration results in the pump speed at any given time 


By specifying the initial pump speed at zero, a pump start transient may be simulated. 

A pump head-flow curve was provided by the pump manufacturer, Barber-Nichols Engineering 
Co. (Fig. 2.3.1). The curve was fitted with a polynomial function to give the head coefficient (y ) 

as a function of the flow coefficient ( (|>) 

^ = 0.52889 - 1.4956(|i + 47.819«ti* - 485.93<|»’ + 1633.9(|i'* - 1833. 5<t>* (2.3.5) 

The flow coefficient (^is obtained &om test data m tenns of the flow rate (in gpm) and the pump 
speed as 


_gpm__ (2.3.6) 

^ 0.0531N 

The pump head is calculated from the pump speed and head coefficient 
H = 4.507 x10‘‘NV (2.3.7) 
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HEAD 






(2.3.8) 


The pump pressure rise is then obtained as 


APp = 


pH 

144 


The jumped pump model requires the pump design flow rate (Q^) and speed (Np ) in order to 
define the other operating characteristics ( HPj, Ip ) required by the model. 


2.4 Tank Thermal Model 

The tank model is a lumped model consisting of four control volumes (Fig. 2.4.i): (1) ullage, (2) 
tank wall, (3) liquid on the tank wall, and (4) bulk liquid. The thermal model of each control 
volume is described in tiie following. 

2.4.1 ilUags 

The ullage thermal nKxlel applies conservation of mass and energy to determine the ullage 
pressure, temperature and mass (Fig. 2.4.2). From conservation of mass, the change in the ullage 
mass ( Mjj ) is due to aU masses entering and leaving the ullage control volume 

(1) droplet evaporation rate in the ullage (mDu) 

(2) boiling rate of the liquid on the tank wall ( mbw) 

(3) bulk liquid boiling rate ( mLu ), or ullage condensation ( muL ) 

(4) liquid surface condensation (mcoND) 


dM„ 

dt 


mou+ mbw+ muj- muL— moc»«D 


(2.4.1) 


These mass flow rates are defined in Section 2.4.6. The ullage mass is obtained by integrating its 
time rate of change with respect to time 




(2.4.2) 


From conservation of energy, the change in the ullage temperature (Tu) is the result of 


(1) heat transfer to the ullage (q^) 

(2) work done on the ullage (w„) 

(3) energy added to the ullage by incoming and leaving masses (ENTHu) 
(j'p Qu ~ ENTHu “ ^vu*^u 


dt 


MyCvU 


(2.4.3) 



m,w *«»!!, 


TANK WALL (W) 


UL I m 


m III m 


L LU 


SPRAY BAR (S) 


Figure 2.4. 1 Tank Thennal Model 
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The terms in Eq. 2.4.3 are defined as follows 

(1) Qu = Qwu “ Quwl - Qul “ Quo " ^us (heat transfer to ullage) 
where 

is the heat-tranfer rate between the tank wall and ullage, 

IqwuI > 0 for a dry wall 
= 0 for a wet wall 

qu^ is the heat-transfer rate between the ullage and wall liquid, 

huwL I = 0 for a dry wall 
> 0 for a wet wall 

q^ is the heat-transfer rate between the ullage and bulk liquid 

q^D is the heat-transfer rate between the ullage and liquid droplet 

q^js is the heat-transferrate between the ullage and (unsubmerged) spray bars 

The above heat-transfer rates are defined in Section 2.4.5. 


(2) wu = p„^i- 


(work done on ullage) 
dV 

where the change in the ullage volume (— ^) is equal and opposite to the change in the liquid 


and wall liquid volumes 


dt 


dt 

(3) ENTH 


dVi. dVwL 


dt 


dt 


(2.4.4) 






where h^ = h^(pu) is the saturated vapor enthalpy of the ullage. 

The ullage volume is obtained as the difference between the tank volume and the bulk liquid and 
wall liquid volumes 

Vl““V^ (2.4.5) 

Eq. 2.4.3 is integrated with respect to time to obtain the ullage ten^)erature 

Tu = (T„),e+/(^)t (2.4.6) 

With the ullage mass, temperature and volume determined, the ullage pressure is calculated from 
the equation of state 


(y- 

i:r: 


£ 


f 

i'- 
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(2.4,7) 


Pu-— 7^ 

2.4.2 Tank Wall 

The tank wall is divided into two sections, one facing the liquid and the other facing the ullage. 
The tank wall facing the bulk liquid is assumed to be at the same temperature as the liquid. Thus, 
the tank wall thermal noodel described in this section applies to the section facing the ullage (Fig. 
2.4.3). Since liquid can form on the tank wall as a result of spraying, the model must account for 
both dry and wet wall cases. 

From conservation of energy, the change in the tank wall temperature is due to 

(1) heat input to the wall finom the environment ( q^w ) 

(2) heat-transfer rate between the wall and ullage ( q^u ) 

|qwu| > 0 for a dry wall 
= 0 for a wet wall 

(3) heat-transfer rate between the wall and liquid on the wall ( q^L)* 

IqwlI = 0 for a dry wall 
> 0 for a wet wall 


— Qew Qwu Qwl (2 4 8) 

dt M^Cpw 

Section 2.4.5 defines these heat-transfer rates. Eq. 2.4.8 can be integrated with respect to time to 
obtain the tank wall temperature 

T.=(T.L+j(^) (2.4.9) 

2.4.3 Wall Liquid 

The wall liquid thermal model is also governed by the laws of conservation of mass and energy 
(Fig. 2.4.4). From conservation of mass, the change in the wall liquid mass ( M^) is ^uai to the 
difference between the liquid mass reaching the wall and the liquid mass boiled off fiom tiie wall 

= msw- mbw (2.4. 10) 

dt 

where msw is the spray flow rate reaching the wall and mbw is the liquid boil-off rate from the 
wall. 

These mass flow rates will be defined in Section 2.4.6. Eq. 2.4.10 can be integrated to obtain the 
wall liquid mass 
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From conservation of energy, the change in the wall liquid temperature ( T^) is the result of heat 
transfer to the wall liquid and sensible energy added to the spray to raise its temperature ( Ts^ ) to 
the wall liquid temperature. Heat transfer to the wall liquid includes heat-transfer rate between the 
wall and wall liquid (q^L)* heat-transfer rate between the ullage and wall liquid (qu^L). 


tiTyyL _ Qwl ~^Quwl °^swCpl(T^ (2 4 12) 

^WL^pWL 

These heat-transfer rates are defined in Section 2.4.5. Eq. 2.4. 12 can be integrated to obtain the 
wall liquid temperature 

Twl = (Tm.).c + 1(^)« (2.4.13) 

The wall liquid vapor pressure is then obtained from the thermodynamic data base as 


PwL ■” pMt(*^WL) 


(2.4.14) 


The volume rate of change of the wall liquid is determined from Eq. 1.9 as 

dV^_ 1 dM^ 
tit dt 

where ^ tiie wall liquid density. 

Eq. 2.4.15 is integrated to obtain the wall liquid volume 

VwL = (^wl)ic J ^ ^ jtit 


(2.4.15) 


(2.4.16) 


2.4.4 Bulk Liquid 

Originally conceived as multi-node, the bulk liquid thermal model is made single node since (1) 
mixing destrat^ the liquid and create a uniform bulk, and (2) uncertainty in heat-transfer 

modeling does not justify the added complexities of a multinode model. 

The liquid thermal model is also based on the laws of conservation of mass and energy. From 
conservation of mass, the change in the liquid mass must be balanced by a change in the ullage 
mass and any mass vented overboard (Fig. 2.4.5). 

dM, 

— : — = msL+ msuL+ mcoND+ muL— muj— ms— mv 
dt 


(2.4.17) 



where msL is the liquid spray flow rate into the bulk liquid 
msuL is the unevaporated droplet flow rate 
mcoND is the liquid surface condensation flow rate 
muLis the ullage condensation flow rate 
iriLu is the liquid bod-off rate 
ms the pump flow rate 
mv is the overboard venting flow rate. 

The liquid mass is obtained by integrating its time rate of change 

= (2.4.18) 

From conservation of energy, the change in the liquid temperature is caused by 

(1) heat transfer to the liquid 

(2) heat added by the unevaporated droplets 

(3) sensible energy added to the liquid spray to raise its temperature ( T$) to the liquid 
temperature 

(4) latent heat of vaporization of the liquid 


(ITl _ + ®SUL Cpi.(Tj “Tl)“ iiiLu(h,j)j^ - liisu Cpl(Tl -Tj) 

4t ^L^pL 

The heat-transfer rate to the liquid ( qj,) is given by 


(2.4.19) 


“ QeL ^UL QlS 

where q^^ is the heat added to the liquid by the environment 

q^j^ is the heat-transfer rate between the ullage and liquid 

qLs is the heat-transferrate between the liquid and (submerged) spray bars 

These heat-transfer rates are given in Section 2.4.5. Eq. 2.4.19 is integrated with respect to time 
to give the liquid temperature 

TL=(Tt),c + j(^)it (2.4.20) 

The liquid vapor pressure is obtained from the thermodynamic data base as 

Pl = P«(Tl) (2.4.21) 

The liquid volume rate of change is determined from the rate of change of the liquid mass 
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(2.4.22) 


dV,. _ 1 dMi. 
dt ~Pl dt 


where Pl = P„,(Tl) is the liquid density. 

Eq. 2.4.22 is integrated to give the liquid volume 

VL=(V0.e+j(^)t (2.4.23) 

2.4.5 Heat Transfer 

This section defines the heat-transfer rates which are found in the energy balances of Section 2.4.1 
to 2.4.4. These heat-transfer rates can be divided into two groups: free convection and forced 
convection. Free convection is the dominant heat-transfer mode in the ullage and liquid, while 
forced convection characterizes liquid droplet heat transfer in the ullage. 

The convection heat-transfer rate is generally defined as 


q = hAAT 


where h is the convection heat-transfer coefficient 
A is the surface area of heat transfer 

ATis the temperature difference between the heat source and sink 
The heat-transfer coefficient is obtained from the Nusselt Number (Nu) as 


h = 




Nu 


where kp is the fluid thermal conductivity and is the surface characteristic length. 


The Nusselt number is a function of the Rayleigh number (Ra) defined as 

aPATL,yc, 

pk 

where a is the acceleration 

P is the tfaennal expansion coefficient, 

P = for gas, -f for liquid 

Tf PVtiT^p 

is the characteristic length 
p is the density 

Cp is the specific heat at constant pressure 
p is the dynamic viscosity 
k is the thermal conductivity 


(2.4.24) 


' 

k' 
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All properties must be evaluated at the film temperature ( Tf ) which is defined as the average of the 
fluid and surface temperatures. 

2.4.5. 1 Free Convection 

Two free convection heat-transfer correlations are used in the model. The first one is a free 
convection correlation for interior surfaces of vertical ducts, vertical plates and cylinders, and 
horizontal cylinders (Ref. 28) (Fig 2.4.6) 

Nu = 0.555Ra®-^ + 0.447 (2.4.25) 

This correlation is used to calculate the heat-transfer coefficients 

( 1 ) between the ullage and wall ( h^ ) 

(2) between the ullage and bulk liquid ( h^i^) 

(3) between the ullage and wall liquid ( hy^ ) 

(4) between the ullage and (unsubmerged) spray bars ( hyj) 

(5) between the bulk liquid and (submerged) spray bars ( hLs) 

The characteristic length for hu^ , h^^and h^^is the internal tank diameter while that of hyg and 
hj^ is the spray bar diameter. 

The second correlation is the McAdams correlation for free convection of vertical surfaces in the 
turbulent range (Ref. 17) 

Nu = 0.13Ra''^ (2.4.26) 

This correlation is used to calculate the heat-transfer coefficient between the wall and wall liquid 
( h^y)- Because of the 1/3 power in Ra, can be obtained without knowing the characteristic 
length, thereby removing the uncertainty in determining the wall liquid layer. 

2.4.5.2 Foregd Convggtioa 

The forced convection heat-transfer coefficient between the ullage and liquid droplets ( hm,) is 
based on a McAdams recommended correlation for flow over a sphere (Ref. 8) (Rg. 2.4.7) 

Nu = 0.3125Re®*“ (2.4.27) 

The Reynolds number of the spray flow (Re) is defined as 

Re = £X^E5a. (2.4.28) 

where Velyis the droplet velocity in the ullage 

Dq is the droplet diameter assumed to be equal to the orifice diameter 
p,pare the density and viscosity of the ullage gas 

Since the droplet diameter and velocity vary with the orifice size, the droplet heat-transfer 
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Rgure 2.4.6 Free Convention Heat-Transfer Correlation for Interior Surfaces of Vertical Ducts, 
Vertical Plates and Cylinders, and Horizontal Cylinders 
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h*D/k 



coefficient must be deteimined for each orifice. The total droplet heat-transfer rate is obtained by 
summing the droplet heat-transfer rates from each orifice 


n 

= S("d-p)i(^d«p)i (2.4.29) 

i=l 

where (nj^). is the number of droplets sprayed fix>m orifice i into the ullage. This is given by 


msu I Dchar 
(n ) - \ A 

V 2pD(VD).(VelD). 

where ^nisu j is the spray flow rate into the ullage from orifice i 

(Vd)j and(VelD) . are the droplet volume and velocity from orifice i 
Pd is the droplet density 

^cHAR is a characteristic length determined empirically. 


(2.4.30) 


By correlating the zero-g TVS model with LeRC ullage pressure collapse data, this characteristic 
length was determined to be 1/4 of the tank diameter. 


2.4.6 Mass Transfer 


This section defines the mass-transfer rates found in the mass balance equations of Section 2.4. 1 to 
2.4.4 which include 

(1) Bulk liquid boiling (muj) 

(2) Liquid boiling from the tank wall ( mt>w) 

(3) Liquid droplet evaporation in the ullage (mou) 

(4) Liquid spray falling into the bulk liquid ( msuL) or acc umula ting on the tank wall ( m$w) 

(5) Ullage condensation (muL> 

(6) Liquid surface condensation (mcoNu) 

2.4.6. 1 Bulk Liquid Boiling 

Bulk liquid boiling occurs when the liqmd vapor pressure is equal to the tank ullage pressure. It 
can be the result of heat transfer to the liquid and/or pressure decay in the ullage. It must also 
include sensible energy added to the liquid spray to increase its temperature to the liquid 
temperature. 

^ Pl ~ Pu» 


mLu = 



— msL c. 


.(T1.-T3)], 


dPu 

dt 


<0 
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(2.4.31) 


{K\ L 


qL-msLCj^(TL-Ts)-MLCpL 


fai^ rdPu 

dt 


^>0 

dt 


A polynomial fit of the LH 2 saturation temperature vs. pressure curve was obtained and its 

9pJ». 


derivative taken to give an expression for 


= 0.37781 -4.9170x10"’pl +21. 7623x10'^ Pl^ 


(2.4.32) 


If the ullage pressure increases above the liquid vapor pressure, boiling stops 
mLu = 0, if Pl<Pu- 
2.4.6.2 Wall Liquid Boiling 

Wall liquid boiling from the tank wall follows the same mechanism as bulk liquid boiling 
^ PwL “ Pu» 


” 7Z T” ^uwL ®sw Cpl(T^ 


)]■ 


^<0 

dt 


where 


K.U 


f— 


qwL+quwL-nwwCpL{TwL-Tsw)-Mwi,CpL^^j . ^>0 


(2.4.33) 


I9p 


= 0.37781 - 4.9170x10'’ p*i + 21.7623x10'^ p,^ 


2 


(2.4.34) 




^ PwL ^ Pu» — 0. 

As with bulk boiling, wall liquid boiling includes heat transfer to the wall liquid and sensible 
energy added to the spray liquid to bring its temperature to the wall liquid temperature. 

2.4.6.3 Liquid Droplet Evaporation in the Ullage 


Liquid droplets in the ullage will start boiling once the subcooled liquid spray is bpught to 
saturation. From an energy balance on the liquid droplets, an expression for the liquid droplet 
boiling is obtained 


mou = 




^UD C 


“ ^s)j 


(2.4.35) 


where Tu^ = T^(pu) is the ullage saturation temperature. 


n-29 



2.4.6A Liquid Soiav Falling into th e Bulk Liquid or Accumulating: on the Tank Wall 

The unevaporated sprayed mass in the ullage is assumed to fall into the bulk liquid under 1 g, or to 
accumulate on the tank wall under 0 g (Fig. 2.4.8), i.e., 


msuL = msu- mou (for 1 g) (2.4.36) 

msw = msu— mou (for 0 g) 

2.4.6.S Ullage Condensation 

Ullage condensation occurs whenever the ullage temperature is equal to the saturation temperature 
corresponding to the ullage pressure. It is the result of heat removal from the liquid droplet (when 
there is spraying) and the wall liquid (Figure 2.4.9) 

T„=T„.(Pu) (2.4.37) 


_ Qud ~^^ul 


2.4.6.6 Liquid Surface Condensation 

When helium is not present to act as a barrier to mass transfer, bulk liquid mixing during pump 
operation induces condensation on the liquid surface. This condensation rate is controlled by the 
heat transfer rate from the ullage to the liquid 


mcoND = 


^UL 

K«)u 


(2.4.38) 
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BULK LIQUID 

Figure 2.4.8 Droplet Evaporation Model 



Figure 2,4.9 Ullage Condensation Model 
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2.6 Sample Cases 

The sample cases shown are for a tank with 90% and 25% liquid quantities and a 0.25 Btu/hr-ft^ 
heat flux. The tank is the 639-ft3 Multi-purpose Hydrogen Test Bed (MTHB) ta^ which is a 
cylindrical tank with elliptical bulkheads at both ends.The cylinder measures 5 ft in length and 10 ft 
in diameter while the bulkhead has a height of 2.5 ft. The tank has a wall thickness of 0.5 in and is 
made of aluminum. One-g acceleration is assumed and no helium is present in the tai^ The 
results show the ullage and liquid vapor pressures, recirculation and vent flow rates, time between 
destratification and venting, destradfication time, and TVS operation fiiequency. 
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Figure 2.6.2 TVS Recirculation Pun^ and Vent Valve Flow Rate Transient During Ullage 

Destratij5cation ( 90 % Liquid Quantity) 
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Figure 2.6.3 TVS Performance Simulation at 25% Liquid Quantity 
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Figure 2.6.4 TVS Performance Simulation at 25% Liquid Quantity 
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SECTION 3 


COMPUTER MODEL DESCRIPTION 


3.1 Programmi ng Description 

The Zero-g TVS performance prediction program was developed on the following system 


Computer 
Operating System 
Language 
Plotter 

Plotting Software 


HP-9000 Series 500 
HP-UX rel. 5.2.1 
FORTRAN 77 rel. 5.12 
HP-7550 
CRTPLT 


3.2 Flow Charts 

3.2.1 Heat Exchan ger Model 

The flow chart of the heat exchanger model is shown in Section 1.4 of Reference 29. 
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3.2.2 Spray Manifold/Iniection Tube Model 
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3.2.3 Integrated Z ero-g TVS Model 
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CALCULATE SPRAY INJECTION FLOW RATE 



CALCULATE HEAT-TRANSFER RATES (WALL/ULLAGE, 
ULLAGEAVALL LIQUID, WALL/WALL LIQUID, 
ULLAGE/BULK LIQUID, ULLAGE/DROPLET,TANK WALL) 


CALCULATE BOIL-OFF RATES (DROPLET, WALL LIQUID, BULK LIQUID), 
AND CONDENSATION RATES (ULLAGE AND LIQUID SURFACE) 


CALCULATE RATES OF CHANGE OF ULLAGE, WALL LIQUID, 
AND BULK LIQUID MASSES 


I 


CALCULATE RATES OF CHANGE OF ULLAGE, BULK LIQUID, 
AND WALL LIQUID VOLUMES 


1 


CALCULATE RATES OF CHANGE OF ULLAGE, BULK LIQUID, WALL 
LIQUID, AND TANK WALL TEMPERATURES 


T 


TIMES TIME + DTIME 


0 



NO 
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3.3 Dftfinition of Variables 

3.3.1 Input Variables 

3.3.1 .1 Heat Exchanger Model 

The input variables of the heat exchanger model are described in Section 3.2.1 of Reference 29. 

3.3. 1.2 Spray Manifold/Iniection Tu be Model 



UNIT 

DESCRIPTION 

dsm 

in 

spray manifold diameter 

zsm 

in 

spray manifold length 

dsi 

in 

spray injection tube diameter 

zsi 

in 

spray injection tube length 

norf 


number of orifices 

nbar 


number of spray bars 

ks 


spray orifice loss coefficient 

cds 


spray orifice discharge coefficient 

rovcrd 


bend r/d of spray manifold 

dmdot 

Ibm/sec 

flow rate increment 

tol 

psia 

convergence tolerance on tank pressure 

nlim 


max number of iterations 

nsec 


number of sections with the same orifice sizes 

node 


node number 

dorf 

in 

orifice diameter 

3.3.L3 Recirculation Pump Model 


SYMBOL 

mn 

DESCRIPTION 
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mdotd 

Ibm/sec 

design pump flow rate 

dpd 

psid 

design pump pressure rise 

npumpi 

rpm 

initial pump speed 

npumpd 

rpm 

design pump speed 

xhp 


multiplier to determine design input horsepower 

xn 


fraction of design pump speed used to determine the 
pump speed operating band 

deltat 

sec 

time needed to reach design speed 

effp 


pump efficiency 

3.3.1.4 Integrated Zero-g TVS Model 


SYMBOL 

uhh 

DESCRIPTION 

xd 


multiplier for spray injection orifice size (to determine 
droplet size) 

xchar 

in 

characteristic length used in the equation to calculate 
the number of droplets 

he 


helium injection indicator (1 = yes, 0 = no) 

xcond 


multiplier for liquid surface condensation rate 

prtsp 


time indicator to print output of subroutine spray 
(output is printed prtsp<time<prtsp+0. 1) 

pui 

psia 

initial ullage pressure 

tui 

R 

initial ullage temperature 

pU 

psia 

initial ullage temperature 

twi 

R 

initial wall temperature 

twli 

R 

initial wall liquid temperature 

fuU 

% 

percent full level 
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xll 

in 

length of straight section downstream 
of recirculation line 90-degree bend 

dtank 

in 

tank diameter 

hcyl 

in 

cylinder height 

hbulk 

in 

tank bulkhead height 

tkw 

in 

tank wall thickness 

dsb 

in 

spray bar diameter 

dl 

in 

diameter of recirculation line upstream of reducer 

62 

in 

diameter of recirculaton line downstream of reducer 

mdvent 

Ibm/sec 

overboard venting flow rate 

mdsi 

Ibm/sec 

initial spray (pump) flow rate 

dthex 

R 

heat exchanger temperature drop 

qflux 

Btu/hr-ft^ 

heat flux 

g 

ft/sec^ 

acceleration level 

pmin 

psia 

control band min pressure 

pmax 

psia 

control band max pressure 

delt2 

sec 

integration time step when pump is off 

ipmt2 


number of time steps between output printing when 
pump is off 

iplot2 


number of time steps between output plotting when 
pump is off 

fintim 

sec 

end time 

deltl 

sec 

integration time step when pump is on 

ipmtl 


number of time steps between output printing when 
pump is on 
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iplotl 


number of time steps between output plotting when 
pump is on 

nline 


number of output lines printed per page 

ovariable 


option to plot variable (l=yes) 

subhd 


plot subheading 

xtitl 


plot x-title 

yvariable 


of variable plot 

3.3. 1.5 LH^ saturation properties 


nsat 


number of data points 

tsat 

R 

saturation temperature 

psat 

psia 

saturation pressure 

enthf 

Btu/lbm 

saturated liquid enthalpy 

shpf 

Btu/lbin~R 

saturated liquid specific heat 

densf 

lWft3 

saturated liquid density 

texpf 

R-1 

saturated liquid thermal expansition coefficient 

condf 

Btu/hr-ft-R 

saturated liquid thermal conductivity 

viscf 

Ibm/ft-sec 

saturated liquid dynamic viscosity 

enthg 

Btu/lbm 

saturated vapor enthalpy 

3.3. 1.6 GHi ProDcrdes 


np (nt) 

number of pressures (temperatures) 

tnrm 


normalized temperature 

tconst 

R 

reference temperature 

pvap 

psia 

pressure 

tvap 

R 

temperature 
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enth 

Btu/lbm 

enthalpy 

shv 

Btu/lbm-R 

specific heat at constant vapor 

shp 

Btu//lbm-R 

specific heat at constant pressure 

dens 

lbm/ft3 

density 

cond 

Btu/hr-ft-R 

thermal conductivity 

vise 

Ibm/ft-sec 

dynamic viscosity 


3.3.2 Output Variables 
3.3.2.1 Heat Exchanger Model 

The output variables of the heat exchanger model are described in Section 3.2.2 of Reference 
29. 


3.3.2.2 Spray Manif old/Iniection Tube Model 


SYMBOL 

mil 

mdot 

Ibm/sec 

dppump 

psid 

dpsm 

psi 

dpsi 

psi 

pin 

psia 

pout 

psia 

pnode 

psia 

ptank 

psia 

mdin 

Ibm/sec 

mdout 

Ibro/sec 

nxls 

Ibm/sec 


DESCRIPTION 

spray (pump) flow rate 

pump pressure rise 

spray manifold tube pressure drop 

spray injection tube pressure drop 

nodal inlet pressure 

nodal outlet pressure 

nodal pressure 

tank pressure 

nodal inlet mass flow rate 

nodal outlet mass flow rate 

spray flow rate through each orifice 
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veld 


ft/sec 


as 

3.3.2.3 Tank Model 

in^ 

SYMBOL 

mn 

pu 

psia 

tu 

R 

vu 


mu 

Ibm 

tw . 

R 

pi 

psia 

tl 

R 

vl 


ml 

Ibm 

pwl 

psia 

twl 

R 

vwl 


mwl 

Ibm 

mv 

Ibm 

ts 

R 

mdlu 

Ibm/sec 

mds 

Ibm/sec 

mdsl 

Ibm/sec 

mdsu 

Ibm/sec 


droplet velocity 
spray orifice area 

PESCRIPTIQN 

ullage pressure 

ullage temperature 

ullage volume 

ullage mass 

wall temperature 

bulk liquid pressure 

bulk liquid temperature 

bulk liquid volume 

bulk liquid mass 

wall liquid pressure 

wall liquid temperature 

wall liquid volume 

wall liquid mass 

mass vented overboard 

spray temperature 

bulk liquid boil-off rate 

spray (pump) flow rate 

spray flow rate into bulk liquid 

spray flow rate into ullage 



mddu 

Ibm/sec 

dropplet evaporation rate 

mdbw 

Ibm/sec 

wall liquid boil-off rate 

mdul 

Ibm/sec 

ullage condensation rate 

mdcond 

Ibm/sec 

liquid surface condensation rate 

qwu 

Btu/sec 

heat-transfer rate between wall and ullage 

quwl 

Btu/sec 

heat-transfer rate between ullage and wall liquid 

qwl 

Btu/sec 

heat-transfer rate between wall and wall liquid 

qui 

Btu/sec 

heat-transfer rate between ullage and bulk liquid 

qud. 

Btu/sec 

heat-transfer rate between ullage and droplet 

qus 

Btu/sec 

heat-transfer rate between ullage and (unsubmerged) 
spray bar 

qls 

Btu/sec 

heat-transfer rate between bulk liquid and (submerged) 
spray bar 

npump 

rpm 

pump speed 

dppump 

psid 

pump pressure rise 
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3.4 Program Listing 

3.4.1 Heat Exchan^erModel 
* 

* Zero G Venting System Integrated Steady State Heat Exchanger 

* Performance Program 

* 

* By Tibor Lak, David Soo Hoo, & Dr. Han Nguyen 

* 

* A HP-9000 Program adapted from the Shuttle Venting Program, Rev. 3 

* Includes single & two phase flow heat transfer & pressure losses. 

* 

* May 8, 1992 

* 

* INITIAL CXDNDrrrrONS: 

« 

* MASSIC - MASS IN THE TANK OR MANIFOLD TO BE VENTED (LBM) 

* PIC -SATURATED PRESSURE OF THE TANK OR MANIFOLD (PSIA) 

* PSTI -INITIAL GUESS AT THE INLET PRESSURE DURING BOILING (PSIA) 

* MDOT -INITALGUESSATTHEVENTORLEAKFLOWRATE(LB/SEC) 

* FLUID PROPERTIES: 

* 

* Prop - Type of Propellant (1.0 = Hydrogen, 2.0 = Oxygen) 

* PTP -TRIPLE POINT PRESSURE (PSIA) 

* AVUG- SONIC VELOCITY IN LIQUID (FT/SEQ 

* 

* EXTERNAL CONDmTONS: 

* 

* PAMB - AMBIENT PRESSURE (PSIA) 

* G -GRAVITY 

* 

* INTERNAL CONDITIONS & CONHGURATIONS: 

* 

* M - NUMBER OF NODES IN THE VENT PATH (Maxium is 20 nodes) 

* VT -TANK OR MANIFOLD VOLUME (CUBIC IN) 

* AX -EXIT AREA (SQ IN) 

* A(l-M) - FLOW AREA OF THE VENT PATH PER NODE (SQ IN) 

* K(l-M)- FLOW LOSS COEFFICIENT PER NODE 

* DH(1-M)- CHANGE IN HEIGHT BETWEEN NODES (IN) 

* QDOT -HEAT FLUX INTO THE FLUID PER NODE (BTU/SEQ 

* (1-M) 

* SAREA - SURFACE AREA OF THE NODE: USED TO CALCULATE THE 

* (1-M) HEAT TRANSFER CX)EFnCIENT 



LENGTH -LENGTH OF NODES (IN) 

* (1-M) 

* 

* VARIABLEEXrr AREA VARIABLES: 

* 

* VA - Change in area (1 = yes , 0 = no) 

* G1 - Acceleration in number of gravities. 

* DIAM - Line diameter (in) 

* SMAX - Maximum exposed surface area (sq inches) 

* 

* PROGRAM CONTROL VARIABLES: 

* 

* FINTIM - MAXIMUM RUN TIME (SEC) 

* PRDEL - PRINT INTERVAL (SEC) 

* QDTERR - QDOT ITERATION ERROR 

* Delptp - Delta Exit Pressure From the Triple Point Pressure 

* DELT -TIME INTERVAL BETWEEN ITERATIONS (SEC) 

* OT - option to plot time vs various parameters (0 = no plot) 

* OM - option to plot mass vs various parameters (0 = no plot) 

* OP - option to plot pressure vs various parameters (0 = no plot) 

* DEBUG - option to print different information between iterations 

* [0.0 = no debug] [ 1.0 = results from mdot loop] 

* [2.0 = results from PST loop] [3.0 = results from PES loop] 

* [4.0 = results from PSFO loop] [5.0 = results from PSFI loop] 

* [6.0 = results from liq PSFO loop] [7.0 = results from PXS loop] 

* [8.0 = results from gas calc, loops] 

* [9.0 = results from all 1 thru 7 loops] 

* 

double precision pic,psti,ptp,avliq,pamb,vt,crrmx,perrmx 
double precision delt,delpq>,ys,sic,time^ot,tsp,tts 
double precision pst,pstl,slt,sgUlrolg:hogt,yUhotv,ptsc,tsl 
double precision berr,ps,pliq,tliq,his,hgs,dmdot,pl,ph 
double precision rhole,ihogc,hlc,hgc,sie,sge,he,qdterr 
double precision rholfo 4 ’hogfo,hlfo,hgfo,slfo,sgfo,hfo»ritofo,pdfo 
double precision htot,pfo,psfo,visgfo,vislfo,asvo,drhdtl,to,qerr 
double precision rholfi^ogfi 4 tifi»hgfi,slfi,sgfi,hfi 4 iiofi,pdfi 
double precision xf,ps^visgfi,vislfi,asvUgxf,phisqf4pff,psfri 
double precision dptf^mf,peir,del,pxs,hlx,hgx,sbc,sgx^obc,yx 
double precision hx,vx jhox^hoxv,pxsc,gerr,tx,pxd,px,thnist 
double precision isp,dpx,sfo 1 Jilf,slf ,psxl /m4p4^hogx 
double precision tdim,ti,tw,pback 

♦ 

double precision pes(2),tes(2),htote(2),pde(2),rhoe(2),hs(2) 

« 

double precision RHOLI(45^),RHOLO(45;2),RHOGI(45,2),RHOGO(45^) 



double precision PSO(45,2),PDI(45,2),PDO(45.2),HI(45,2),HO(45,2), 
1 TSI(45,2),tave(45,2) 

double precision RHOO(45,2),SI(45,2),SO(45,2),pi(45,2),po(45,2), 

1 tso(45,2) 

double precision VI(45,2),VO(45,2),MACHI(45,2),MACHO(45,2) 
double precision DPF(45^),PHISQ(45,2),DPT(45,2),DPM(45,2), 

1 ERR(45,2),rhoi(45,2) 

double precision XTT(45,2),QDTT(45,2),RN(45,2),psi(45^) 
double precision Xl(50),p2(50),p3(50),p4(50),p5(50),p6(50),p7(15) 
double precision p8(50),p9(50),p 10(50), lmxpar(50)4visp(50) 
double precision lrhop(50),grhop(50),rhovp(50),lent(50),lentp(50) 
double precision gent(50),sclent(15,15),pl l(50),pl2(50),drhdto(50) 
double precision t7(15,15),gvisp(50),gentp(50) 
double precision cpg(6,14),thkg(6,14),tl4(14),viscg(6,14), 

1 pl4(6),delta(45) 

realyo(45,2),yi(45,2),a(45,2)4c(45,2),qdot(45,2),dh(45,2) 
reaimass 4 rf,niassic,mlow 4 nhigh,mdot( 2 ),twall( 45 ),ye( 2 ),mdMn 8 
character *8 name( 10 ),subt( 10 ) 
integer exit,ot,om,op,num(15),va 

DeHne Heat Transfer Variables 

real length(45),sarea(45,2) 

double precision pl3(50),surft(50),q(45,2),cpl(50),thkl(50), 

1 stemp(50),cp,ts,st,cond 

Define Plot Variables 

real ptitle(18,4) 

real subhd(l 8 ) 4 iame 1 (1 0 ) 4 iame 2 ( 10 ),name 3 (l 0),name4(10) 
real name 6 ( 10 ) 4 iamc 7 ( 10 ),name 8 ( 10 ) 4 iame 9 ( 10 ),namel 0 ( 10 ) 
real namell( 10 ) 4 iaincl 2 ( 10 ),name 5 ( 10 ) 
conimon /pltcom/ niisc(3),nc4niss(13),dclirn,ltick,nfig,nptnrtin, 

1 nlines,nchlin,ptitle 

data name/'mdot ','isp ’,’thrust *,'qerr ’ 

1 ,'rhot ','pst ’/pcs ',’htot ' 

2 ,'ext area'/twali 7 

Read Input Data 
read (5,100) lable 

read (5,101) massic,pic,psti,mdot(l),pback 

read (5,100) lable 

read (5,101) prop,ptp,avliq,pamb,g 



read (5,100) lable 

read (5,102) m,vt,ax,va,tw 

read (5,100) lable 

read (5,*) mdot(2),disp,dth,vtdi,ed 

read (5,100) lable 

read (5,101) dpinc,errmx,perrmx,debug,exdi 
read (5,100) lable 

read (5,101) fintim,prdel,qdterr,delt,delptp 

read (5,100) lable 

read (5,104) ot,om,op,dq 

read (5,100) lable 

do i=l,m 

read (5,105) a(i,l),dh(i,l),qdot(i,l),length(i) 
enddo 

read (5,106) subhd 
read (5,107) subt 

100 format (//,al) 

101 format (5fl0.0) 

102 format (il 0 , 2 fl 0 . 04 l 0 ,fl 0 . 0 ) 

103 format (4fl0.0,el0.1) 

104 format (3il0,fl0.0) 

105 format (4fl0.0) 

106 format (/,18a4) 

107 format (10a8) 

* 

* Write Input Data 

if (debug .eq. 1.0 .or. debug .ge. 10.0) then 
write (6,101) massic,pic,psti 4 iidot(l),pback 
write (6,101) prop,ptp,avliq,pamb,g 
write (6,102) m,vt,ax,va,tw 
write (6,*) mdot(2),disp,dth,vtdi,ed 
write (6,101) dpinc,eirmx,pemnx,debug,exdi 
write (6,101) fintim,prdel,qdterr,delt,delptp 
write (6,104) ot,om,op,dq 
do i=l,m 

write (6,105) a(i,l),dh(i,l),qdot(i,l)4ength(i) 
enddo 
endif 

♦ 

* read tables of data file misc.data 
if (prop .eq. 1.0) then 

open (unit=2, file*’h2misc.data',status= unknown’) 
else 
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open (unit=2, file='o2niisc.data',status='unknown') 
endif 

read (2,150) nl 
do i=l,nl 

read (2,155) xl(i),lmxpar(i) 
enddo 

* 

read (2,150) n2 
do i=l,n2 

read (2,155) p2(i)4visp(i) 
enddo 

* 

read (2,150) n3 
do i=l,n3 

read (2,155) p3(i),gvisp(i) 
enddo 

* 

if (prop .ne. 1.0) then 
read (2,150) nl2 
do isl,nl2 

read (2,155) pl2(i),drhdto(i) 
enddo 
endif 

* 

* read tables of data file iho.data 
if (prop .eq. 1.0) then 

open (unit=3, file*'h2rho.data',status= unknown’) 
else 

open (unit=3, file='o2iho.data',status='unknown') 
endif 

read (3,150) n4 
do i=l,n4 

read (3,155) p4(i)4rhop(i) 
enddo 

* 

read (3,150) n5 
do i=l,n5 

read (3,155) p5(i),griiop(i) 
enddo 

4c 

read (3,150) n6 
do i=l,n6 

read (3,155) riiovp(i),p6(i) 
enddo 
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read (3,151) n7,maxt 
do i=l,n7 

read (3,*) p7(i),num(i) 
doj=l,num(i) 

read (3,155) t7(i,j),sclent(i,j) 
enddo 
enddo 

read tables of data file enthalpy 8c entropy data 
if (prop.eq.1.0) then 

open (unit=s4, file=’h 2 entdata‘,status=’unknown’) 
else 

open (unit=4, file='o2entdata’,status='unknown') 
endif 

read (4,150) n8 
do i=l,n8 

read (4,155) p8(i),lentp(i) 
enddo 

read (4,150) n9 
do i=l,n9 

read (4,155) p9(i),gentp(i) 
enddo 

read (4,150) nlO 
do i=l,nl0 

read (4,155) plO(i),lent(i) 
enddo 

read (4,150) nil 
doi-l,nll 

read (4,155) pll(i),gcnt(i) 
enddo 

read tables of data file thermal conductivity & surface tension data 
if (prop.eq.1.0) then 

open (unit=10, file='h2thermo.data’,status= unknown') 

read (10,150) nl3 
do i=l,nl3 

read (10,*) pl3(i),cpl(i),thkl(i),suift(i),stemp(i) 
enddo 



read (10,154) nl4,numl,tcnst 
do i=l,nl4 
read (10,152) pl4(i) 
doj=l,numl 

read (10,153) cpg(i,j),thkg(ij),viscg(i,j) 
enddo 
enddo 

doj=l,numl 
read (10,152) tl4(j) 
enddo 
endif 

150 format (//43) 

151 format (//,i3,4x43) 

152 format (flO.2) 

153 format (fl4.0,2el4.3) 

154 format ^/,i3,4x43,4x/10.0) 

155 format (2f 14.0) 

* 

* Define Program Constants & working variables 

* 

ys = 0.0 
AI«A(M,1) 

AE = A(1,1) 
tsp = PTP 
K(M,1)»0.5 
K(l,l) = 0.5 
PES(l) = FTP 
PES(2)*PTP 
ml = m - 1 
mass ~ massic 

call value(nl0,pl04ent,pic,sic) 

beta » 1.0 

gamma si.o 

time = 0.0 

m8 * m * 0.8 

* pltime =s -pltdel 
prdme » -prdel 
ncountl = 1 
ncount2 s 1 
exit =s 0 

pyi = 3.14159 

* 

do ij * 1, m 
twall(ij) = tw 




a(u,2) = pyi * disp ** 2.0 / 4.0 
sarea(ij,l) = pyi * (disp + 2.0 * dth) * length(ij) 
sarea(ij^) = pyi * disp * length(ij) 
qdot(ij,2) = -qdot(ij,l) 
enddo 

Define Plot Constants 


nc = 0 
nfig = 0 
nlines = 0 

Check to Determine the Correct Time increment is Used 

(deltgtpltdel) delt = pltdel 
(deltgt.prdel) delt * prdel 

CALCULATE INITIAL TANK OR MANIFOLD CONDITIONS 

1RH0T = MASSAT’ 
tagl *0.0 
tagh = 0.0 
mlow = 0.0 
mhigh = 0.0 
pst = PSTI 
pstl = 0.0 
jl=0 


Loop to Calculate Initial Pst 

5 call valuc(nl0,pl04cntJPST,slt) 
call value(nl l,pl l,gent,PST,sgt) 
call value(n4,p44rhopJPST4iiolt) 
call value(n5,p5,grhop4*ST^hogt) 

YT « (SIC-SLT)/(SGT-SLT) 
if (yt .It 0.0) yt * 0.0 
if (yt .gt 1.0) yt = 1.0 

RHOTV = YT/(1.0/RHOT-1.0/RHOLT*(1.0-YT)) 
call value(n64iovp,p6,RHOTV,ptsc) 
if (jl .eq. 0)pstl =ptsc 
beta = dabs(pst • pstl) 
pstl s pst 
beir « ptsc - pst 
IF (dabs(beir) .gt 0.001) then 
if (berr .gt 0.0) then 


pst = pst + beta/2.0 
else 

pst = pst - beta/2.0 
endif 

if (jl.gt.lOO) go to 6 
go to 5 
endif 
6 PS=PST 

IF (debug .eq. 2.0 .or. debug .ge. 10.0) THEN 
write (6,110) jl,pst,slt,sgt4‘holt,rhogt,ptsc,mass,yt 
1 10 format ('counter = ',i4,4x,'PST = ',f8.3y,7(3x48.3)) 
ENDDF 

PSTI=PST 
PLIQ=PST 
if (prop.eq.1.0) then 
call H2SAT(PST,diq) 
else 

call o2sat(PST,tliq) 
endif 
YS=YT 


SATURATED UQUID AT SOURCE 

IF(PLIQ.le.PST) then 
call value(n8,p84entp,PSTJils) 
call value(n9,p9,gentpJPST,hgs) 

HS(1) = (1.0-YS)*HLS+YS*HGS 
HS(2) = HS(1) 
else 

SUBCOOLED LIQUID AT SOURCE 
YS=0.0 

call value3(np7,maxMium,p7,t7,sclent,pliq,tliq,hls) 
HS(1) = HLS 
HS(2) = HS(1) 
endif 


Set Constants for Mdot Calculation 



c START ITERATION ON THE FLOWRATE CONVERGENCE BASED ON CHOKE 
FLOW 1 

c AT THE ORinCE BY THE OUTLET OF THE TANK FOR THE TVS SYSTEM 


THETA = 1.05 
DO LI = 1, 1000 

if (tagl .ne. 0.0 .and. tagh .ne. 0.0) then 
mdot(l) = mlow + dm * (pi / (pi - ph)) 
else 

mdot(l) = (theta + 1.0) ♦ mdot(l)/2.0 
endif 

* definition OF TANK EXIT CONDITION (SONIC/SUBSONIC FLOW) 

* (CHOKE FLOW AT THE TANK EXIT NODE M) 

HTOTE(1)=HS(1)+QDOT(M,1)/MDOT(1)-G*DH(M,1)/(12.0*1728.0) 

AT = A(M,1) 
scold = 0.0 
xpes = 0.9 
tsp - xpes * tsp 
do i = 1,400 
tsp ss tsp + dpinc 
IF(tsp.GT.PST) GO TO 15 
IF(tsp.LT.PTP) tsp = ptp 
call valuc(n4,p4,lriiop,tsp,rhole) 
call valuc(n5,p5,grhop,tsp,rhoge) 
call value(n8,p84entp,tsp,hle) 
call value(n9,p9,gentp,tsp,hge) 
call value(nl0,pl0,lent,tsp,slc) 
call value(nl l,pl l,gcnt,tsp,sgc) 

caUTPHS(HTOTE(l),RHOLE,RHOGE,MDOT(l)AT.HLE,HGE,SLE, 
1 sge,ye(l),he,se,rhoc(l),pde(l)) 

if (se .Ic. seold) goto 15 
seold = se 
enddo 

15 continue 

if (prop.cq.1.0) then 
caUH2SAT(tsp.tts) 
else 

call o2sat(tsp,tts) 
endif 


* 


PO(M,l)=tsp+PDE(l) 

YO(M,l)=YE(l) 



PRESSURE LOSS CALCULATION THROUGH NODE M 


KF=K(MJ) 

AREA=A(M,1) 

CHANGE IN ELEVATION 

DHI=DH(M,1) 

DHO=DH(M+l,l) 

HO(M,1)=HS(1)+QDOT(M,1)/MDOT(1)-G*DH(M4)/(12.0*1728.0) 

HTOT=HO(M.l) 

TWO-PHASE FLOW REGION 

PFO=PO(M,l) 

PSFO=PFO 

DETERMINE THE OUTLET CONDITIONS AT NODE N, GIVEN THE TOTAL 
PRESSURE AT THE OUTLET 

DO II = 1.5 

IF(PSFO.LT.PTP) PSFO=PTP 
IF(PSFO.GT.PST) PSFO=PST 
call valuc(n4,p4,lihop,PSFO.rholfo) 
call value(n5,p5,grhop,PSFO.ihogfo) 
call value(n8,p84entp,PSFO,hIfo) 
call value(n9,p9,gentp,PSFO,hgfo) 
call value(nlO,plO,lent,PSFO,slfo) 
call value(nll,pll,gent.PSFO,sgfo) 
call valuc(n3,p3,gvisp,PSFO,visgfo) 
call valuc(n2,p24visp J*SFO,vislfo) 

caUTPHS(HTOT4UIOLFO,RHOGFO,MDOT(l),AREAJILFO.HGFO,SLFO, 
1 sgfo,yfo4ifo.sfo.rhofo.pdfo) 

PSFO=PFO-PDFO 

End of II Loop (To Determine Node M Outlet Conditions) 
enddo 

if (prop.eq.1.0) then 
call H2SAT(PSfo,tsl) 
else 

call o2sat(PSfo.tsl) 
endif 


PSO(M,l)=PSFO 

PDO(M,l)=PDFO 

RH0L0(M,1)=RH0LF0 

RHOGO(M,1)=RHCXjFO 

RH00(M.1)=RH0F0 

YO(M,l)=YFO 

H0(M,1)=HF0 

SO(M,l)=SFO 

VO(M,1)=144.0*MDOT(1)/(A(M,1)*RHCXXM,1)) 

TSO(M,l)=tsl 
if (prop.le.1.0) then 

caU H2SVEL(YFO,PSFO,RHOLFO,RHOGFO,asvo) 
else 

call value(nl2,pl2,drhdto,psfo4rhdtl) 
call o2sat(psfo,to) 

call o2svel(yfo,to4>olfojhogfo,drhdtl,asvo) 
endif 

MACHO(M,l)=VO(M,l)/ASVO 
IF(MACHO(M,1).LT.O.O) MACHO<M,1)=0.0 
IF(MACHO(M,1).GT.1.0) MACHO(M,1)=10 

INITIAL GUESS AT NODE INLET BASED ON COMPRESSIBLE LOSS 

call INITL(PSFO,KF,MDOT(l),AREA,RHOLFO,psa) 

TOTAL ENTHALPY AT NODE INLET 

HI(M,1)=HS(1)4QDOT(M+1,1)/MDOT(1)-G»DH(M.I)/(12.0*1728.0) 

HTOT=HI(M,l) 

DEFINITION OF STATIC PRESSURE AT NODE INLET 
DOI2 = 140 

IFO*SFLLT.PTP) PSFI=PTP 

IF(G.LE.O.OJ^.PSFLLT.PSFO) PSFI=PSFO 

IF(PSFLGT.PST) PSFI=PST 

call value(n4,p4 Jihop,PSFl4holfi) 

call value(n5,p5,griiopJ*SFI^hogfi) 

call value(n8,p8Jentp^SFIJilfi) 

call value(n9,p9,gentp,PSFI,hgfi) 

call value(nl0,pl04ent,PSFI,slfi) 

call value(nll,pl l.gent^SFI^gfi) 

call value(n3,p3,gvisp,PSFI,visgfi) 

call value(n2,p2Jvisp JSFI,vislfi) 

call TPHS(HTOT4IHOLFIJIHOGnj4DOT(l)AREA3LFI,HGFI,SLFI, 
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sgfi,yfi,hfi,sfi4’hofi,pdfi) 

PFI=PSFI+PDFI 

* LOCKHART-MARTINELU TWO-PHASE FLOW PARAMETER 

* 

caUXPARAM(YFO,VISLFO,VISGFO,RHOLFO,RHOGFO,Yn,VISLn, 

1 visgfi,RHOLH,RHOGn,xf) 

IF(XF.LT.0.01) XF=0,01 

LGXF=LOG10(XF) 

call value(nl,xl,lmxpar,LGXF,phisqf) 

PfflSQF=(10.0**PfflSQF)**2 

♦ 

* TWO-PHASE FLOW PRESSURE LOSS AND NODE INLET PRESSURE (PSFN) 

* 

caUTPS(AREA,G,DHI,DHOJCF,PSFO,YFO,RHOLFO,RHOGFO,YFI, 

1 rholfi,RHOGFI,MDOT(l)^HISQF,dpff,psfn,dptf, 

2 dpmf) 

ERR(M,1)=ABS(PSFN-PSFI) 

PSFI=(PSFI+PSFN)/2.0 

* 

* CONDITION WHERE THE INLET STATIC PRESSURE OF THE NODE IS 

DETERMINED 

♦ 

IF(ERR(M,1).LE.ERRMX) go to 20 

* 

* End of 12 Loop (To Determine Node N Inlet Conditions) 


enddo 

20 PFI*PSFI+PDn 

* 

* CALCULATE THE INLET CONDITIONS OF THE NODE M 

* 

RHOU(M,l)=RHOLFI 

RH0GI(M,1)=RH0GFI 

RHOI(M,l)=RHOFI 

DPM(M,1)*DPMF 

DPF(M,1)*DPFF 

DPT(M4)=DPTF 

PI(M,1)=PH 

PDI(M,1)=PDFI 

PSI(M,1)=PSFI 

YI(M,l)=Yn 

HI(M,1)=HH 

si(M,i)=sn 

PHISQ(M,1)=PHISQF 



VI(M,1)=144.0*MDOT(1)/(A(M,1)*RHOI(M,1)) 
if (prop.le.1.0) then 

call H2SVEL(YH,PSFI,RH0LFI,RH0GFI,asvi) 
else 

call value(nl2,p 12,drhdto,psfi,drhdt2) 

call o2sat(psfi,ti) 

call o2svel(yfi,d^holfi^hogfi,drhdt2,asvi) 
endif 

MACHI(M,1)=VI(M,1)/ASVI 
IF(MACHI(M,1).LT.0.0) MACHI(M,1)=0.0 
IF(MACHI(NU).GT.1.0) MACHI(M,1)=1.0 

CALCULATE THE TWO PHASE HEAT TRANSFER COEFFICIENT FOR NODE M 
BASED ON INLET PRESSURE 

call value(nl3,pl3,cpl,psff,cp) 
call value(nl3,pl3,thkl,psfi,cond) 
call value(nl3,pl3,surft,psfi,st) 
call value(nl3,pl3,stemp,psfi,ts) 
if (tave(m,l) .le. 0.0) then 
tavc(m,l) = ts 
endif 

hgl = hgfi - hlfi 
dtemp = twall(m) - tavc(m,l) 
di^exdi 

dprcss = 778.2 ♦ dtemp * rhogfi * hgl / ts 
if (yfi .gt 0.7) then 

re = 48.0 * mdot(l) * (1.0-0.7)/(pyi * vislfi * di) 
else 

re = 48.0 * mdot(l) ♦ (1.0-yfi)/(pyi * vislfi * di) 
endif 

callhtcoeff(phase4temp4pmss,cp»cond,visM^holfi, 

1 rhogfi,st,hgUc^,(ii,qdott) 

if (yfi .ge. 0.7) then 
tdim = 0.0 

call value2(nl4,numl,pl4,tl4,cpg,psfi,tdim,cp) 
call value2(n 14^uml ,p 1 4,tl4,thkg,psfi,tdim,cond) 
call value2(nl4,numl ,p 14,tl4,viscg,psfi,tdim,visgfi) 
re = 48.0 * mdot(l) /(pyi * visgfi * i) 
call htcoeff(3.0, dtemp, dpress,cp,cond,visgfi^holfi, 

1 rhogfi, st,hgl4©,xf,di,qtt) 



qdott = qdott*(l-yfi) + qtt*yfi 
endif 

q(m4) = qdott * (twall(m) - tave(m,l)) * sarea(m,l)/144.0 

TSI(M,1)=TS 

XTT(M,1)=XF 

RN(M,1) = RE 

QDTT(M,1) = QDOTT 

TAVE(M,1) = (tsi(m,l) + tso(m,l))/2.0 

theta = ps/pi(m,l) 
perr = ps-pi(m,l) 
if (perr .It 0.0) then 
del = mlow - mdot(l) 
tagl = 1.0 
mlow = mdot(l) 
pi = perr 
else 

del = mhigh - mdot(l) 
tagh = 1.0 
mhigh = mdot(l) 
ph = perr 
endif 

dm = mhigh - mlow 

if (dabs(perr) .le. perrmx .or. dabs(del) .le. 0.0001) goto 30 
enddo 


c END OF MOOT CALCULATION LOOP I 


30 continue 
c 

c Determine the Back Pressure to the Heat Exchanger part of the Spray 
c Bar 
c 

ptk = psti 
mdt =: mdot(2) 
dspry ^disp 

call spray(mdt,ptk,dspry,ed,pback) 

♦ START ITERATION ON THE HEAT TRANSFER COEFFICffiNT ♦ 

DOL=1,150 
if (1 .gt 1) then 
do ii = l»ml 
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nn = m-ii 

qdot(nn,l) = qdot(nn+l,l) + q(nn,l) 
qdot(nn,2) = qdot(nn+l,2) - q(nn,2) 
if (sarea(u,l) .le. 0.0 .or. qdtt(ii,l) .It 0.0) then 
cl = 1.0 
else 

cl = sarea(ii,2) * qdtt(ii,2)/(sarea(ii,l)*qdtt(ii,l)) 
endif 

if (dabs(delta(ii)) .gt dq) then 
tw = twall(ii) 

twall(ii) = (cl * tave(u,2) + tave(ii,l))/(l+cl) 
twall(ii) = (tw + twall(ii))/2.0 
if(twall(U) .It 24,845) twall(ii) = 24.845 
en^ 
enddo 
endif 
qdt = 0.0 
qdtl = 0.0 

c START THE VENT & SPRAY NODAL NETWORK FLOW PROPETIES 
CALCULATION I 


do il = 1, 2 

* 

* DEFINITION OF EXIT CONDITION (SONIC/SUBSONIC FLOW) 

* (CHOKE FLOW AT THE EXIT NODE) 

HKyrE(U)=HS(U)+QDOT(l,a)/MDOT(U)+G*(DH(M+141)-DH(l4l))/ 
1 (12.0*1728.0) 

AE = A(141) 
scold = 0.0 
xpes = 0.9 

PES(il) = xpes*PES(U) 
do i = 1^00 
if (il .cq. 1) then 
PES(il) = PES(U) + dpinc 
else 

PESCil)*pback 

endif 

IF(PES(il).GT.PST) GO TO 40 
IF(PES(U).LT.PTP) PES(il)=PTP 
call valuc(n4,p4dihop,PES(il)^ole) 
call value(n5,p5,grhop,PES(il),rhoge) 
call valuc(n8,p8,lentp,PES(il),hle) 
call value(n9,p9,gentp,PES(il),hge) 
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call value(nl0,pl04ent,PES(il),sle) 
call vaiue(nl l,pl l,gent,PES(il),sge) 

callTPHS(HTOTE(il),RHOLE,RHOGE,MDOT(il),AE,HLE,HGE,SLE, 

1 sge,ye(il),he,se 4 'hoe(il),pde(il)) 

if (se .le. scold) goto 40 
scold = sc 
cnddo 

40 continue 

* 

IF (debug .eq. 3.0 .or. debug .ge. 10.0) THEN 
write (6,111) i,pes(il),mdot(il),ye(il),pde(il),sle,sge, 

1 rhole,rhoge,se,seold 

111 format Ccounter = *44,4x,'PES = ',f8.3,4x,'MDOT = 

1 fl0.5,2(f8.3,3x)y,6(3x,f8.3)) 

ENDIF 

if (prop.eq. 1 .0) then 
call H2SAT(PES(il),tes(il)) 
else 

call o2sat(PES(il),tes(il)) 
endif 

* 

♦ 

PO(l,il)=PES(il)+PDE(U) 

YO(l,il)=YE(U) 

C+4.+++++.. H .++-h H- f I I I I H f - M f f H -++++++++++- H i -+- H -++4 f f + +++-M I M M I M f H+ + + 
++++ 

c PRESSURE LOSS CALCULATION THROUGH EACH NODE + 

c (NODEIIS AT THE EXIT & NODE MIS AT THE INLET) + 

1 f M ^ f f +f4++4 H -++++-K i H -f +++++4- ++++ 

++++ 

DON=l,M 

if(il .eq. 1 .and. n .eq. m) goto 60 
IF(N.GT.l) PO(N41)=PI(N-141) 

IF(N.GT.l) YO(N41)=YI(N-141) 

IF(N.GT.l) TSO(N41)=TSI(N-l,il) 

AREA=A(N41) 

« 

if (N .eq. 1 .and. PSFO .Ic. 0.0) psfo = po(n41) 
call value(n2,p24visp,PSFO,vislfo) 
vise = vislfo 
if (il .eq. 1) then 
if (n .ne. 1) then 
dd = vtdi/12.0 
call Mct(dd,cd,mdot(l),visc,ff) 

K(n,l) = ff * length(n)/vtdi 
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else 

K(U) = 0.5 
endif 
else 

dd = disp / 12.0 

call frict(dd,ed,mdot(2),visc,ff) 

K(n,2) = ff * length(n)/disp 
endif 

KF=K(N,U) 

* 

* CHANGE IN ELEVATION 

♦ 

DHI=DH(N,U) 

DHO=DH(N+141) 

HO(N41)=HS(U)+QDOT(N,U)/MDOT(il)+G*(DH(M+141)-DH(N,U))/ 
1 (12.0*1728.0) 

HTOT=HO(N,U) 

IF(YO(N,il).LE.O.O) PHASE^l.O 
IF(YO(N,il).GT.O.O) PHASE=2.0 

* 

* DETERMINE THE PHASE OF THE FLUID 

* 

IF(PHASE.ge.2.0) then 


* TWO-PHASE FLOW REGION 


PFO=PO(N,U) 

PSFO=PFO 

* 

* DETERMINE THE OUTLET CONDITIONS AT NODE N, GIVEN THE TOTAL 

* PRESSUREATTHEOUTLET 

* 

DO II = 1,5 

IF(PSFO.LT.PTP) PSFO=PTP 
IF(PSFO.GTJPST) PSFO=PST 
call valuc(n 4 ,p 44 rhop,PSF 04 ’holfo) 
call value(n5,p5,grhop,PSFO,riiogfo) 
call value(n8,p84en^ J^SFO,hlfo) 
call value(n9,p9,gentp J*SF04igfo) 
call value(nl0,pl04ent,PSFO,slfo) 
call value(nll,pll,gent,PSFO,sgfo) 
call value(n3,p3,gvisp J*SFO,visgfo) 
call value(n2,p24visp J*SFO,vislfo) 

caU TPHS(HTOTJlHOLFO,RHOGFO,MDOT(il),AREA,HLFOjaGFO, 

1 slfo,sgfo,yfo,hfo,sfo,rhofo,pdfo) 
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PSFO=PFO-PDFO 


End of II Loop (To Determine Node N Outlet Conditions) 
enddo 

IF (debug .eq. 4.0 .or. debug .ge. 10.0) THEN 
write (6,112) n,il,psfo,slfo,sgfo,rhofo,yfo,sfo,hfo,pdfo 
1 12 format (’node #'42,4x,'counter = ',i4,4x,'PSFO = 

1 f8.3y,7(3x,f8.3)) 

ENDEF 

PSO(N,U)=PSFO 
PDO(N41)=PDFO 
RHOLO(N,il)=RHOLFO 
RHOGO(N,U)=RHOGFO . 

RHOO(N41)=RHOFO 

YO(N,U)=YFO 

HO(N41)=HFO 

SO(N,il)=SFO 

VO(N41)«144.0*MDOT(il)/(A(N,il)*RHOO(N4l)) 

calculate the outlet temperature at node 1 

if (n .eq. 1) then 
callH2SAT(psfo,tsl) 
call value(nl3,pl3,cpl,psfo,cp) 
hgl = hgfo - hlfo 
qq = qdot(141)/mdot(il) 
if (yo(141) .gt 0.99) then 
if (il .eq. 1) then 

tso(l,l) = tsi(m,l) + (qq - hgl)/cp 
else 

tso(l,2) = tliq + (qq + hgl)/cp 
endif 
else 

tso(141) -tsl 
endif 
endif 

if (prop.le.1.0) then 

call H2S VEUYFO J>SFO,RHOLFO,RHOGFO,asvo) 
else 

call value(nl2,pl2,drhdto,psfo,drhdtl) 
call o2sat(psfo,to) 


call o2svcl(yfo,to/holfo^hogfo,drhdtl ,asvo) 

endif 

MACHO(N,U)=VO(N,il)/ASVO 
IF(MACHO(N,il).LT.O.O) MACHO(N,il)=0.0 
IF(MACHO(N,U).GT.1.0) MACHO(N,U)=1.0 

★ 

* INITIAL GUESS AT NODE INLET BASED ON COMPRESSIBLE LOSS 

* 

call INlTL(PSFO JCF>lDOT(il),AREA,RHOLFO,psfi) 

* 

• TOTAL ENTHALPY AT NODE INLET 

4c 

HI(N41)=HS(a)+QDOT(N+141)/MDOT(U)+ 

1 G*(DH(M,U>DH(N+141))/(12.0*1728.0) 

HTOT=HI(N,U) 

4c , - 

* DEFINITION OF STATIC PRESSURE AT NODE INLET 

4c 

DO 12 = 1,50 

IF(PSFLLT.PTP) PSn=PTP 

IF(G.LE.O.O,AND.PSFI.LT.PSFO) PSFI=PSFO 

IF(PSn.GT.PST) PSn=PST 

call value(n4,p4,lriiop,PSFI,rholfi) 

call value(n5,p5,grhop,PSFI,rhogfi) 

call value(n8,p84entp,PSFI,hlfi) 

call value(n9,p9,gentp4*SFI,hgfi) 

call value(nlO,plO,lent,PSFI,slfi) 

call value(nl l,pl l,gent,PSFI,sgfi) 

call value(n3,p3,gvisp4*SFI,visgfi) 

call value(n2,p24visp4>SFI,vislfi) 

call TPHS(lfrOT4RHOLFI41HOGn,MDOT(U),AREA,HLFI,HGFI, 

1 slfi,sgfi,yfi4'fi*sfi,rliofi*pdfi) 

PFI=PSFI+PDFI 

4c 

♦ LOCKHART-MARTINELU TWO-PHASE FLOW PARAMETER 

4c 

callXPARAM(YFO,VISLFO,VISGFO,RHOLFO,RHOGFO.YFI,VISLFL 
1 visgfi,RHOLFLRHOGFLxf) 

IF(XF.LT.0.01) XF=0.01 

LGXF=LOG10(XF) 

call value(nlpcl4mxpariGXF,phisqf) 

PHISQF=(10.0**PHISQF)*»2 

4c 

* TWO-PHASE FLOW PRESSURE LOSS AND NODE INLET PRESSURE (PSFN) 
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caU TPS(AREA,G,DHI,DHO,KF,PSFO,YFO,RHOLFO,RHOGFO,YH, 

1 rholfi,RHOGFI^DOT(il),PHISQF,dpff,psfn,dptf, 

2 dpmf) 

ERR(N,il)=:ABS(PSFN-PSFI) 

PSn=(PSFI+PSFN)/2.0 

* 

* CONDITION WHERE THE INLET STATIC PRESSURE OF THE NODE IS 

DETERMINED 

* 

IF(ERR(N,U).LE.ERRMX) GO TO 50 

* 

* End of 12 Loop (To Determine Node N Inlet Conditions) 

♦ 

enddo 

50 Pn=PSFI+PDH 

* 

IF (debug .eq. 5,0 .or. debug .ge. 10.0) THEN 
write (6,113) n42,psfi,slfi,sgfi,rhofi,yfi,sfi,hfi, 

1 phsiqf 

113 format (*node #‘42,4x,'counter = ’,i4,4x,'PSFI = 

1 f8.3y,6(3x,f8.3)y,3x,f8.3) 

ENDIF 

* CALCULATE THE INLET CONDITIONS OF THE NODE N 

♦ 

RHOLI(N,U)=RHOLH 

RHOGI(N41)=RHOGFI 

RHOI(N41)=RHOFI 

DPM(Nal)=DPMF 

DPF(N41)=DPFF 

DPT(N41)=DPTF 

PI(N,il)=PFI 

PDI(N,il)=PDFI 

PSI(N,U)=PSn 

YI(N41)=YFI 

HI(N41)=Hn 

SI(N,U)=Sn 

PHISQ(N,ii)=PHISQF 

VI(N41)=144.0*MDOT(il)/(A(N,il)*RHOI(N,il)) 
if (prop.le.1.0) then 

call H2SVEL(YnjPSFI,RHOLFI,RHOGFI,asvi) 
else 

call value(nl2,pl2,drhdto,psfi,drhdt2) 
call o2sat0psfi,ti) 

call o2svel(yfi,ti,rholf!,rhogfi,drhdt2,asvi) 



endif 

MACHI(N,il)=VI(N.il)/ASVI 
IF(MACHI(N,a).LT.O.O) MACHI(N,U)=0.0 
IF(MACHI(N,a).GT.1.0) MACHI(N,a)=1.0 

CALCULATE THE TWO PHASE HEAT TRANSFER COEFFICIENT FOR NODE N 
BASED ON INLET PRESSURE 

call value(nl3,pl3,cpl,psfi,cp) 
call value(nl3,pl3,thkl,psfi,cond) 
call value(nl3,pl3,surft,psfi,st) 
call value(nl3,pl3,stemp,psfi,ts) 

do nnl = 1, 2 

if (tave(n41) .le. 0.0) then 
tave(n41) = ts 

endif 

hgl = hgfi - hlfi 

if (il .eq. 1) then 
dtemp = twall(n) - tave(n41) 
di~exdi 
else 

dtemp = tave(n,il) - twall(n) 
di = disp 
endif 

dpress = 778.2 * dtemp * rhogfi * hgl / ts 
if (yfi .gt 0.7) then 

re = 48.0 * mdot(U) * (1.0-0.7)/(pyi * vislfi * di) 
else 

re = 48.0 ♦ mdotOl) ♦ (1.0-yfi)/(pyi * vislfi * di) 
endif 

write (6,*) dtemp,dpicss,cp,cond,vislfi;rholfi^hogfi, 

1 st4igl,re^ 

caUhtcocff(phase,dtemp4prcss,cp,cond,vislfi4tolfi, 

1 rhogfi,st,hgl^,xf,di,qdott) 

if (yfi .ge. 0.7) then 
tdim = 0.0 

callvalue2(nl4,numl,pl4,tl4,cpg,psfi,tdim,cp) 
call value2(n 144iuml ,p 14,tl4,thkg,psfi,tdim,cond) 
call value2(nl44iuml ,p 14,tl4,viscg,psfi,tdim,visgfi) 
re = 48.0 ♦ mdot(il) /(pyi ♦ visgfi * di) 



1 


callhtcoeff(3.0,dtemp,dpress,cp,cond,visgfi^holfi, 
rhogfi,st,hgl^,xf,di,qtt) 
qdott = qdott*(l-yfi) + qtt*yfi 
endif 

if (il .eq. 1) then 

q(n,il) = qdott * (twall(n) - tave(n41)) * sarea(n,il) 

1 /144.0 

else 

q(n,il) = qdott ♦ (tave(n41) - twall(n)) * sarea(n,il) 

1 /144.0 

endif 

calculate the inlet temperature at node il 

call value(nl3,pl3»cpl,psfi,cp) 
if (yi(n,il) .gt 0.99) then 
if (il .eq. 1) then 

tsi(n41) = tso(n 41 ) - q(n,il)/(mdot(il) ♦ cp) 
else 

tsi(n41) = tso(n41) + q(n,il)/(mdot(il) * cp) 
endif 
else 

tsi(n,il) s ts 

endif 

XTT(N,il) = XF 
RN(N,U) = RE 
QDTT(N,il) = QDOTT 
TAVE(N41) = (tsi(n,il) + tso(n,il))/2.0 

enddo 

if (debug .ge. 10.0) dien 
write (6,120) tsi(n 41 )«twall(n),tave(n 41 ) 4 %,qdott, 
1 q(n,il),saiea(n,il) 

120 format(4(2x,fl4.4)y,3(2x,fl4.4)) 

endif 
ELSE 


SINGLE PHASE FLOW (PHASE=1.0, LIQUID FLOW) 


CALCULATION OF OUTLET CONDITIONS GIVEN TOTAL PRESSURE AT NODE 
OUTLET (PO(N)) 


PFO=PO(N,il) 



PSFO=PFO 
DO 13 = 1,5 

IF(PSFO.LT.PTP) PSFO=PTP 
IF(PSFO.GT.PST) PSFO=PST 
call value(n4,p4,lrhop,PSFO,rholfo) 
PDFO=144.0*MDOT(U)/A(N,il)*MDOT(il)/A(N,il)/ 
(2.0*RHOLFO*32.2) 

PSFO=PFO-PDFO 

End of 13 Lx)op (To Determine Node N Outlet Conditions) 


enddo 

IF(PSFO.LT.PTP) PSFO=PTP 
IF(PSFO.GT.PST) PSFO=PST 
call value(n8,p8,lentp,PSFO,hfo) 
call value(nlO,pl04entJ*SFO,sfol) 
call value(nl3,pl3,cpl,psfo,cp) 
call value(n2,p24visp4*SFO,vislfo) 


if (n .eq. 1) then 
if (il .eq. 1) then 

tso(141) = tsi(m,il) + qdot(l,il)/(mdot(il) * cp) 
else 

tso(l,il) = tiiq + qdot(l,il)/(mdot(il) * cp) 
endif 
endif 

PSO(N,U)=PSFO 

PDO(N41)=PDFO 

RHOLO(N,il)=RHOLFO 

RHOGO(N,il)=0.0 

RHOO(N,il)=RHOLFO 

YO(N,il)=0.0 

HO(N41)=HS(il)+QDOT(N41)/MDOT(il)+G* 

(DH(M+141)-DH(N,U))/(12.0*1728.0) 

S0(N41)=SF01 

VO(N41)=144.0*MDOTCil)/(A(N,il)*RHOO(N41)) 

ASVO=AVLIQ 

MACHO(N41)«VO(N,U)/ASVO 
IF(MACHO(N41).LT.O.O) MACHO(N,il)=0.0 
IF(MACHO(N41).GT.1.0) MACHO(N,U)=1.0 


IF (debug .eq. 6.0 .or. debug .ge. 10.0) THEN 
write (6,114) n43,psfo,pdfo,tso(n41),riiolfo,hfo, 
so(n,il),vo(n,il),htote(il),tso(n41) 



114 format ('Node # = ',i4,4x,'Counter = ',i4,4x, 

1 'Liquid PSFO = ',f8.3y,8(3x,f8.3)) 

ENDIF 

* 

» INLET CONDITION CALCULA'HONS 
* 

RHOLI(N,il)=RHOLO(N41) 

RHOGI(N,il)=0.0 
RHOI(N41)=RHOU(N41) 

DPM(N,U)=0.0 

DPF(N,a)=144.0*K(N,il)/(2.0*RHOLI(N,il)*32.2)* 

1 (MDOT(U)/A(N,il))**2 

DPT(N41)=DPF(N41) 

PI(N4I)= PO(N,a)+DPF(N 4 l)-G*RHOLI(N,U)» 

1 (DH(N,U)-DH(N+l4l))/1728.0 I 

PDI(N41)=144.0*(MDOT(il)/A(N,il))**2/ -ve 

1 (2.0*RHOU(N,U)*32.2) [I 

PSI(N,il)=PI(N41)-PDI(N41) 

YI(N41)=0.0 

VI(N41)=144.0*MIX)T(il)/(RHOLI(N,U)*(A(N,il))) 

ASVI=AVLIQ 

MACHI(N41)=VI(N,U)/ASVI 

PSFI=PSI(N,il) I- 

IF(PSn.LT.PTP) PSFI=PTP ^ 

IF(PSFI.GT.PST) PSFI=PST 

call value(n2,p24visp,PSFI,vislfi) 

call value(n8,p84entp,PSFI,hlf) 

call value(nl0,pl04ent^Sn.slf) 

* HI(N,U)=HLF 

HI(N41)=HS(il)4QDOT(N+141)/MDOT(U)+ ^ 

1 G*(DH(M41)-DH(N+141))/(12.0*1728.0) 

SI(N,U)=SLF 

PHISQ(N,il)=1.0 

c 

IF (debug .eq. 6.0 .or. debug .ge. 10.0) THEN j 

write (6^14) n,psfi,pdi(n41),dpni(n41), 

1 dpf^(ii41),dpt(n4l),vi(n41) 

214 format CNodc # » *44,4x,’Liquid PSH = ’,f8.3y, > 

1 5(3x,f8.3)) 

ENDIF 

3k i . 

I 

* CALCULATE THE HEAT TRANSFER COEFnCIENT FOR NODE N 

* BASED ON INLET PRESSURE £ 

* 

can value(nl3,pl3,cpl,psfi,cp) 
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call value(nl3,pl3,thkl,psfi,cond) 
call value(nl3,pl3,surft,psfi,st) 
call value(nl3,pi3,stemp,psfi,ts) 
if (tave(n,il) .le. 0.0) tfien 
tave(n,il) = tso(n,il) 

endif 

hgl = hgfi - hlfi 

if (il .eq. 1) then 
dtemp = twall(n) - tave(n,il) 
di = exdi 
else 

dtemp = tave(n,il) - twall(n) 
di = disp 
endif 

dpress « 778.2 * dtemp * rhogfi * hgl / ts 
re = 48.0 * mdot(il)/(pyi * vislfi * di) 

callhtcoeff(phase,dtemp,dprcss,cp,cond,vislfi4iiolfi, 
1 rhogfi,st,hgl^pcf,di,qdott) 

if (il .eq. 1) then 

q(n,il) = qdott * (twall(n) - tave(n41)) * sarea(n,il) 

1 /144.0 

else 

q(n,il) = qdott * (tave(n41) - twaU(n)) * sarea(n41) 

1 /144.0 

endif 

if (t1 .eq. 1) then 

tsi(n41) = tso(n 41 ) - q(n41)/(mdot(il) * cp) 
else 

tsi(n41) = tso(n,il) + q(n 41 )/(mdot(il) * cp) 
endif 

XTT(N,il)=XF 
RN(N,il) = RE 
QDTT(N41)* qdott 
TAVE(N 41) = (tsi(n,il) + tso(n,il))/2.0 

IF (debug .eq. 6.0 .or. debug .ge. 10.0) THEN 
write (6,314) tsi(n,il),q(n41),dtemp,xf,re,qdott, 

1 hgl,cp,di,tave(n,il) 

3 14 format (’liquid TSI =* 'i8.3,3xil0.4,2(3x,f8.3),3x, 
1 fl2.2y,5(3x,f8.3)) 

ENDIF 
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ENDIF 

* PO(N+l,U)=PI(N,il) 

C+4-+++ <-'i - f+4-++++++-< - H+++++++-f-f+-f-f-l-H-+ f ) t i"f-f-+4’+++++- H -+++++ +++++++++ - f H -++ 
+++++ 

c End of N do loop (Conditions for Both the Inlet & Outlet of the + 
c Nodal Network are Defined) + 

C++++++++++4-+++++++-f+++++4 ” f++4-++-f‘+++’f-4'- H -+++-f4 -H -+ + f +-l- + f 4 f ++-h++++-f t + f 4-+ 
++++ 4 - 

60 continue 
ENDDO 

c 

c RE-CALCULATE THE NODAL PROPERTIES OF THE NODES THAT ARE 100% GAS 

c 

if (il .eq. 1) then 
do ig as ml, 1, -1 
if (yo(ig,l) .gt 0.99) goto 130 
enddo 

130 phase = 3.0 

dummy = 1.4 ♦ 766.55232 ♦ 32.2 
do ia = ig, 1, -1 
if (yo(ia,l) .ge. 0.99) then 
iia = ia + 1 
psi(ia,l) = pso(iia,l) 
tsi(ia,l) = tso(iia,l) 
hi(ia,l) = ho(iia,l) 

rhogi(ia,l) =s psi(ia,l) ♦ 144.0 / (766.55232 ♦ tsi(ia,l)) 
rholi(ia,l) = 0.0 
rhoi(ia, 1) = ihogi(ia,l) 

vi(ia,l) = 144.0 * mdot(il) / (a(ia,l) * rhoi(ia,l)) 
asvi = (dummy * dabs(tsi(ia,l)))**0.5 
machi(ia,l) = vi(ia,l)/asvi 
if (debug .eq. 8.0 .or. debug .ge. 10.0) then 
write(6,9990) ia,psi(ia,l),tsi(ia, l),hi(ia, 1 ) 
write(6,9991) vi(ia,l),asvi 4 nachi(ia,l) 
write(6,9992) ihoi(ia,l),rholi(ia,l),ihogi(ia,l) 

9990 formatCnode # =',3x433x,'ps, ts, h (in) =',3(2x,fl0.4)) 

9991 formatCvel, a, m (in) =’3(2x,fl0.4)) 

9992 formatCrho(in) =’,3(2x,fl0.4)) 
endif 

* 

* CALCULATE THE HEAT TRANSFER COEFRCIENT FOR 100% GAS NODES 

* BASED ON INLET PRESSURE 

* 

call value(nl3,pl3,stemp,psi(ia,l),ts) 
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tdim = (tsi(ia,l)-ts)/(tcnst-ts) 
if (tdim .It 0.0) tdim = 0.0 
call value2(nl4,numl,pl4,tl4,cpg,psi(ia,l),tdini,cp) 
call value2(nl4,numl,pl4,tl4,thkg,psi(ia,l),tdim,cond) 
call value2(nl4,numl,pl4,tl4,viscg,psi(ia,l),tdim,visgfi) 
ti = tsi(ia,l) * (1.0 + 0.2 * machi(ia,l) ** 2.0) 
pi(ia,l) = psi(ia,l) ♦ (1.0 + 0.2 * machi(ia,l) ** 2.0) 
pdi(ia,l) = pi(ia,l) - psi(ia,l) 
hgl = hgfi - hlfi 

♦ 

dtemp = 0.0 
dpress = 0.0 

re = 48.0 * mdot(il) /(pyi * visgfi * exdi) 

callhtcoeff(phase,dtemp,dpress,cp,cond,visgfi^holfi, 

1 rhogfi,sthgl^,xf,exdi,qdott) 

if (debug .eq. 8.0 .or. debug .ge. 10.0) then 
write(6,9993) cp,cond,visgfi 
write(6,9994) ti,pi(ia,l),pdi(ia,l)4«,qdott, 

1 twall(ia),tave(ia,l) 

9993 formatCcp, k, mu =’,3(2x^l0.4)) 

9994 formatCTo, Po, Pd(in) =’,3(2x^l2.4)y, 

1 ’Re, He ,Twall. Tave ^',4(2x,fl2.4)) 
endif 

* 

* CALCULATE THE OUTLET PROPERTIES OF THE GAS NODE 

* 

do nnl = 1,2 

if(tave(ia,l) .le. 0.0) tave(ia,l) = ti 

♦ if(nnl .eq. 1) tave(ia,l) = ti 
dtemp = twall(ia) - tavc(ia,l) 

* 

to = ti + (pyi * length(ia) * qdott ♦ (twall(ia)-ti) * 

1 exdi / (mdot(il) ♦ cp * 144.0)) 

if (to Jt 0.0) to = ti/2.0 

call htandf(machi(ia, l),ti,to,twall(ia),macho(ia, 1)) 
ratio = (l+0.2*machi(ia,l)**2.0)/(l+0.2*macho(ia,l)**2.0) 
tso(ia,l) = tsi(ia,l) * (to/ti) * ratio 
pso(ia,l) = psi(ia,l) ♦ (machi(ia,l)/macho(ia,l)) ♦ 

1 (tso(ia,l)/tsi(ia,l))**0.5 

vo(ia,l) = vi(ia,l) ♦ (machi(ia,l)/macho(ia,l)) ♦ 

1 (tso(ia,l)/tsi(ia,l))**0.5 

po(ia,l) = pi(ia,l) * (tso(ia,l)/tsi(ia,l)) * ratio**3.5 
pdo(ia,l) = po(ia,l) - pso(ia,l) 
ho(ia,l) = hi(ia,l) + cp * (to-ti) 
riioo(ia,l) = rhoi(ia,l) * pso(ia,l) ♦ tsi(ia,l) / 
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1 


(psi(ia»l) * tso(ia,l)) 
rholo(ia,l) = 0.0 
rhogo(iaj) =rhoo(ia,l) 
q(ia,l) = qdott * dtemp * sarea(ia,l)/144.0 
tave(ia,l) = (tsi(ia,l) + tso(ia,l))/2.0 
qdtt(ia,l) = qdott 
if (debug .eq. 8.0 .or. debug .ge. 10.0) then 

write(6,9995) pso(ia,l),to,tso(ia,l),ho(ia,l),rhoo(ia,l), 

1 vo(ia,l),asvi,macho(ia,l),q(ia,l), dtemp, 

2 twall(ia),tave(ia, 1 ) 

9995 formatCPs, To, Ts(out) =',3(2x,fl2.4)y,*H, rho, vel(out) =’ 

1 ,3(2x,f 12.4)y,’a, M, Qdot (out) =‘,3(2x,f 12.4)/ 

2 ’Dt , TwaU ,Tave =’,3(2x,fl2.4)) 
endif 

enddo 

endif 

m(ia,l) = re 
enddo 
endif 


c END OF RE-CALCULATING GAS NODES 

c 

enddo 

c 

c END OF THE VENT & SPRAY NODAL NETWORK FLOW PROPETIES 
CALCULATION I 


* 

* CALCULATE THE TOTAL QDOT & NEW TWALL FOR ALL THE ELEMENTS IN 
THE NETWORK 

flagl = 0.0 
dojj = l,m 
qdt = qdt + q(jj,l) 
qdtl = qdtl + q(jjj,2) 
delta(ij) « q(ij,l)-q(ij^) 
if (dabs(delta^» ie. dq) flagl = flagl + 1.0 
enddo 

qcrr = qdtl - qdt 

* qenr = qdot(l,l) - qdt 

* 

IF (debug .eq. 1.0 .or. debug .ge. 10.0) THEN 
write (6,116) l,time,mass,mdot(l),ps,qerr,pes,po(l,l), 

1 po(2,l),po(3,l),po(4,l),po(ll,l),po(12,l), 

2 pi(l,l),pi(2,l),pi(3,l),pi(4,l),pi(ll,l). 


m-40 



3 pi(12,l) 

1 16 format (’iteration # ',i3,2x,'time =',f8.3,3x,'mass =’, 

1 f8.3,4x,'mdot =',fl0.5y,'pst =’,f8.3,4x,'qerr 

2 f8.3,4x,'pes =',f8.3y,6(3x,f8.3),/,6(3x,f8.3)) 

ENDIF 

if (debug .eq. 9.0) then 
write(6, 117) l,qdt,qdtl ,qerr 
dojj = l,m 

write(6,l 18) jj,twall(j|j),q(jj4),q(jj,2),delta(ij), 

1 tsi(ij,l),tso(jj,l),tsi(ij,2),tso(jj,2), 

2 qdtt(ij,l),qdtt(jj^) 
enddo 

1 17 format ('iteration # ’,i3,2x,’Gas Qdot =',f8.3,3x, 

1 ‘Liq Qdot ='^.4,3x,’Q error =’,f8.4y/,t3,’Node #*, 

2 tl3,’T wall',t23;Gas Q',t33,’Liq Q',t43;Delta Q’,t53, 

3 'G Tsi',t63,’G Tso',t73/L Tsi’,t83,'L Tso',t93,'G He’, 

4 tl03,’L HdJ) 

118 format (6x44,10(2x,f8.4)) 
endif 

* CX)NDmON TO END THE QDOT CONVERGENCE LOOP ♦ 

* * 

IF (dabs(qerr) .le. qdterr .or. flagl .ge. m8) goto 200 
* * 

* End of L do loop (The Calculated Qdot to the Nodal ♦ 

* Model is Within the Error Margin of PERRMX * 

* to the past Qdot) ♦ 

* 

ENDDO 

200 do ii s Lml 

nn = m-ii 

qdot(nn,l) = qdot(nn+l,l) + q(nn,l) 
qdot(nn,2) = qdot(nn+l,2) - q(nn,2) 
enddo 
c 

c DETERMINE THE EXIT PLANE CONDITIONS 
c 

j2 = 0 

if (va .eq. 1) ax = sa 
PXS=PSO(l,l) 

205 IF(PXS.LT.PTP) PXS^PTP 
IF(PXS.GT.PST) PXS=PST 
call value(n8,p8dentp,PXS,hlx) 
call value(n9,p9,gentp,PXS,hgx) 
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call value(nlO,plO,lent,PXS,slx) 
call value(nl l,pl l,gcnt,PXS,sgx) 
call value(n4,p4,lrhop,PXS^holx) 
call value(n5,p5,grhop,PXS,rhogx) 
YX=(SO(U)-SLX)/(SGX-SLX) 
HX=YX*HGX+(1.0-YX)*HLX 
VX=SQRT(2.0*32.2*777.649*ABS(HTOTE(1)-HX)) 
RHOX=144.0*MDOT(1)/(VX*AX) 
RHOXV=YX/(1.0/RHOX-1.0/RHOLX*(1.0-YX)) 
call value(n6^hovp,p6,RHOXV,pxsc) 
if (j2 .eq. 0) psxl = pxsc 
gamma = dabs(pxs - psxl) 
psxl = pxs 
gerr = pxsc - pxs 
IF (dabs(geir) .gt .0.001) then 
if (gerr .gt 0.0) then 
pxs = pxs + gamma/2.0 
else 

pxs = pxs - gamma/2.0 
endif 
j2=j2+l 

if (j2.gtl50) go to 210 
go to 205 
endif 

t 

210 IF (debug .eq. 7.0 .or. debug .ge. 10.0) THEN 

write (6,115) j2,pxs,time,gerr,mass,sbc,sgx,rhox,yx,sx,hx 

1 15 format Ccounter = '44,4x,'PXS * ',f8.3,3(3x,f8.3)y, 

1 6(3x,f8.3)) 

ENDIF 

if (prop.eq.1.0) then 
callH2SAT(pxs,tx) 
else 

call o2sat(pxs,tx) 
endif 

PXD=144.0/(2.0*RHOX*32.2)*(MDOT(1)/AX)*(MDOT(1)/AX) 
PX=PXS+PXD 
FM=MDOT(l)*VX/32.2 
FP=PXS*AX 
THRUST=FM+FP 
ISP=THRUST/MDOT(l) 
dpx = pxs - ptp 
flag = 0.0 
flagl = 0.0 
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flag3 = 0.0 

IF (dpx .le. delptp) then 
exit= 1 
ENDIF 

* 

* Print Output 

do 11 = 1, 2 
if (11 .eq. 1) then 
write (6,1000) subhd,subt 
write (6,1430) 

write (6,1100) name(l),mdot(ll),name(2),isp,name(3),thrust, 

1 name(4),qeir,name(5)4*hot,name(6),pst,name(7), 

2 pesGl)»name(8),htot,name(9),a(l,ll) 
write (6,1200) 

write (6,1300) pxs,tx,rhox,yx,hx,vx 
write (6,1400) pes(ll),tesGl)»ihoeGl),yeGl).htote(ll) 
else 

write (6,1000) subhd,subt 
write (6,1460) 

write (6,1100) name(l),mdotGl),name(2),isp,name(3),thnist, 

1 name(4),qerr4iame(5)>rhot,name(6),pst,name(7), 

2 pesGl)»name(8),htot,name(9),a( 1 ,11) 
write (6,1200) 

write (6,1400) pesGl),tesGl),ihoeGl),yeGl),htoteGl) 
endif 

do ncl = l,m 

write (6,1500) ncl,pso(ncl,ll),tsi(ncl,ll),rhoo(ncl,ll), 

1 yo(nc l,ll),ho(nc 1 ,ll),vo(ncl,ll),qdot(nc 1,11), 

2 a(nc m)»phisq(nc 1 41),length(nc l),dpf(nc 1 4 1)» 

3 dpm(ncl41),dpt(ncl,ll),xtt(ncl41)^(ncl,ll), 

4 qdtt(ncl,ll) 
enddo 

write (6,1600) pst,tiiq,iiiolt,yt,hsGl) 
enddo 

write (6,1000) subhd,subt 
write (6,1700) ps,tsp,ttS4ndot(l) 
do 11 = 1, 2 
if G1 -cq. 1) then 
write (6,1430) 
write (6,1800) 
else 

write (6,1460) 
write (6,1800) 
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endif 

do nc2 = l,m 

write (6,1900) nc2,pi(nc2,ll),psi(nc2,ll),pdi(nc2,ll), 

1 po(nc24 1 ),pso(nc2,l 1 ) ,pdo(nc2,l 1 ),tsi(nc2,l 1 ), 

2 tso(nc2,ll),twall(nc2),tave(nc2,ll),q(nc2,ll), 

3 delta(nc2)4c(nc2,ll) 
enddo 

enddo 

ic 

1000 formatClVA18a4y,10a8,/) 

1100 format(9(4x,a8,f9.4y)) 

1200 format(t2,'NODE #*, til, ’STATIC F,t23,'TEMF,t32,’DENSITY’, 

1 t42,‘QUALITY’,t51,’EmTIALPr,t61,’VEL<XITY’,t73,’QDOT, 

2 t83,'AREA’,t93,'PfflSQ’,tl02,’LENGTH’,tll4,’DPF,tl24,'DPM’, 

3 tl34,’DFT,tl44,’Xtt',tl53,’Re #’,tl60,'Ht Trans C’) 

1300 formate/, t2, 'OUTLET, tl 149.4,4(x,f9.4),x,f9.2) 

1400 format(/,t2,’EXrr,tl I,f9.4,3(x,f9.4),x,f9.2) 

1430 format(/,t2,’FLOW PROPERTIES OF THE VENT PART OF THE SYSTEM’) 
1460 format(/,t2,’FLOW PROPERTIES OF THE HEAT EXCHANGER PART OF, 
1 • THE SYSTEM’) 

1500 format(t443,tl l,f9.4.7(xjE9.4)pc,f9.2,5(x,f9.4),x49.2,x49.4) 

1600 formate/, t2,’INLET,tll,f9.4,5(x,f9.4)) 

1700 format(/,t2,’Ptank static', t21,f9.4y,t2,’Pchoke’,t21,f9.4y,t2, 

1 Tchoke’,t21,f9.4y,t2,’Mdot',t2149.4) 

1800 format(/,t2,’NODE #‘,tl3,’Pt in',t25,’Ps in’,t34,’Pd in’, 

1 t44,’Pt out’,t53,’Ps out',t63,'Pd out’,t73,Ts in’,t83, 

2 Ts out’,t93,T waU’,tl03,T ave’,tl 13,’Qdot’,tl23, 

3 ’DelQ’,tl33,*Kloss'y) 

1900 format(t443,tll,f9.4,12(x49.4)) 

STOP 

END 

subroutine H2SAT(PST,tsat) 

% 

double precision pst,tsat,^,t 

« 

FP=LOG(PST/187.506) 

T=L00003+FP*(2.12094E-l+FP*(2.83129E-2+FP*1.75686E-3)) 

TSAT=59.3568*T 

return 

END 

subroutine h2SVEIXY,PS JlSL4lSV,as) 

double precision ps 4 sljrsv,as,p 4 pf»tstl,ts,tsm,tc,tfun,t,tk 
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double precision pa,dpvdt,d2pvtl,d2pvdt,csattl,csatt2,csat 
double precision drodtl,tt,drodt2,drodt,rho,dsdrho,dsdpl,dsdp 
c 

P=PS/187.506264 

FPF=LOG(P) 

TSTl=FPF*(2.83128E-2+FPF*(1.75686E-3)) 

TS=1.00002+FPF*(2. 12094E- 1+TSTl) 

TSM=TS*32.976 

TC=32.976 

TFUN=(TC-TSM)**0.33333 
T=TSM*1.8 
IF(Y.LE.0.0) Y=0.0 
TK=T/1.8 

PA=10.0**(2.00062-50.0970/(TK+1.0044)+1.74849E-2*TK)*14.696 

DPVDT=PA*2.30258*(50.0970/(TK+1.0044)**2+L74849E-2)/1.8 

D2PVT1=DPVDT**2/PA 

D2PVDT=D2PVTl+PA*2.30258/1.8**2*(-2.0*50.0970/(TK+1.00044)**3) 

CSATTl=1.68157*TK/(32.976-TK)**0.1-32.8027 

CSATT2=TK*(3.35743E-2+TK*(-7.68297E-4+6.90292E-6*TK)) 

. CSAT=(CSATT1+TK*(6.81698+TK*(-0.731943+CSATT2)))/2.01572 

DRODT1=-0.38* 7.32346E-3/(32.976-TK)**0.62+4.40742E-4 
TT=(32.976-TK)**0.333333 

DRODT2*TT*TT*(1.66667*2.92263E-04-2.0*4.00849E-05*rr) 

DROLDT=(DRODT1-1.33333*6.62079E-4*TT+DRODT2)*69.9099 

RHO=1.0/(1.0/RSL+Y*(1.0/RSV-1.0/RSL)) 

DSDRHO=-1.0/RHO/RHO*DPVDT*144.0 

DSDP1=777.649/144.0*CSAT/(T*DPVDT)+1.0/RSL/RSL*DROLDT 

DSDP=DSDP1+(1.0/RHO-1.0/RSL)*D2PVDT/DPVDT 

AS=SQRT(ABS(-DSDRHO/DSDP*32.2)) 

return 

END 

subroutine htcoeff(phase,dt,dp,cpJc,mu^hol 4 iiov,st 4 ambda^,x, 

1 diazn,qdot) 

c 

c This subroutine calculates the two phase heat transfer coefficent 
c 

double precision cp4t,mu^ol»rhov,st,x 
double precision Igx(12),lgf(12),lgsf(12),sf(12),f,s,c2,d3 
real lambda 
c 

data lgx/-1.0, -0.824,-0.745,-0.699,-0.585, 0.0,0. 176,0.398,0.663, 

1 0.778,1.778,2.0/ 

datalgf/0.0,0.0212,0.0414,0.0569,0.1139,0.4472,0.5563,0.699, 

1 0.8751,0.9542,1.699,1.8692/ 
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data lgsf/0.176, 0.415, 0.672, 1.0, 1.19, 1,29,1.398, 1.544,1.602, 

1 1.663,1.778,2.0/ 

data sf/0.84,0.73,0.58, 0.37,0.25,0,2,0. 16,0. 12,0. 1 1,0. 1,0.09, 

1 0.09/ 

c 

dtemp = abs(dt) 
dpress = abs(dp) 
gc = 32.2 
al = .023 
a2 = rc**0.8 

a3 = (cp ♦ mu * 3600. / k)**0.4 
a4 = (k* 12.)/(3600. ♦diam) 
hfc = al * a2 ♦ a3 * a4 
if (phase .It 1.5 .or. phase .ge. 3.0) then 
qdot = hfc 

c write(6,*) cp,muJc,diam,al,a2,a3,a4,qdot 
return 
endif 

c write(6,*) al,a2,a3,a4 
c 

bl = .00122 

b2 = ((k/3600.)**0.79) ♦ (cp**0.45) * (rhol**0.49) * (gc**0.25) * 

1 (dtemp**0.24) ♦ (dpress**0.75) 

b3 = ((st*12.)**0.5) ♦ (mu**0.29) * Gambda**0.24) ♦ (rhov**0.24) 
hfz«bl*b2/b3 
c write(6,*) bl,b2,b3 
c 

cl = 1 /x 
c2 = logl0(cl) 
if(c2 .gt.2.0)c2 = 2.0 
call value(124gx4gf,c24) 
f = 10.0**f 
c write(6,*) cl,c2/ 
c 

dl=re*f*1.25 
d2«dl/ 10000.0 
d3 = logl0(d2) 
if(d3.gt.2.0)d3=2.0 
call value(124gsf,sf,d3,s) 
c write(6,*) dl,d2,d3 
qdot = hfc*f + hfz*s 
c write (6,100) b2,b3,cl,c2,dl,d2 

c write (6,100) rc,hfc44tfz,s,qdot 

clOO format (3(3x,fl4.4)y3(3x/14.4)) 
return 
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subroutine INlTL(PO,K,MDOT,A,RHO,pin) 
double precision po,rho,pin4pf 

real k,mdot 

DPF=144.0*K*((MDOT/A)**2.0)/(2.0*RHO*32.2) 

PIN=PO+DPF 

return 

END 

subroutine spray(mdot,psmi,dsni»ed,pback) 
subroutine spray.f to model flow in the spray injection tube 
character*! label 

double precision rod( 12 ),lod( 12 ) 4 'overd,loverd,pback 
real mdin(200),mdout(200),mds(200),as(200),vel(200) 
real pin(200),pout(200),pnode(200),ptank(200),x(200),delpt(200) 
real node(10),asec(10) 

real subhd( 1 8)^titl( 1 8),yp( 1 8),ymdot( 1 8),y mds( 1 8), 

1 yvel(18),yas(18) 

real mdot,mdoti»mdoto 4 ndots,mdotp 
real kfsnU^bsm^ccsm^cbsUcs 
integer*2 iboxdloc 
common /contr]Abox 41 oc 

datanb/12/rod/L,L5,2.,3..4.,6..8.,10.,12.a4.,16.,20./ 
datalod/20.,14.,12.,12.,14..17.,24.,30.,34.,38.,42.,50./ 
ibox = 1 
iloc = 0 


open (unit=8, flles'spray.data*,status='unknown') 

read (8,100) label 

read (8,*) pu,zliq,ztank,dmdot 

read (8,100) 

read (8,*) zsm,dsi,zsi,n,nbar 
read (8,100) 

read (8,*) ks,cds,acc,rho,visc,roverd 
read (8,100) 
read (8,*) tol^ilim 
read (8,100) 

read (8,*) opn,opt,omdin,omds,ovel,oas 
read (8,101) subhd,xtitl,yp,ymdot,ymds,yvel,yas 



read (8,*) nsec 
do i = l,nsec 

read (8,*) node(i),asec(i) 
enddo 

100 format (//a 1) 

101 format (18a4/18a4/18a4/18a4/18a4/18a4/18a4) 
write (6,1) 

1 format (5x,’SPRAY MANIFOLD AND INJECTION TUBE FLOW MODEL’/) 
write (6,2) zliq,pu,psmi,acc,rfio,zsm,dsm 

2 format (5x,'Liquid Level = ',f6. 1 ,' in'/ 

1 5x,’UUage Pressure = ’,f6.2,’ psia'/ 

2 5x,'Spray Manifold Inlet Pressure = ',f6.3,' psia’/ 

3 5x,’ Acceleration Level = ’,f6.1,’ g’/ 

4 5x,’Liquid Density = ’,f6.3,’ lbm/ft3’/ 

5 5x,’Spray Manifold Tube Length » ’,f6. 1,’ in’/ 

6 5x,’Spray Manifold Tube ID = ’4^6.2,’ in’) 
dzsi = zsi/n 

write (6,3) nbar,zsi,dsi,n,dzsijcs 

3 format (5x,’Number of Spray Injection Tubes =* ’,i6/ 

1 5x,’Spray Injection Tu^ Length = ’,f6.1,’ in’/ 

2 5x, ’Spray Injection Tube ID = ’,f6.3,’ in’/ 

3 5x,’Number of Orifices = ’,i6/ 

4 5x, ’Orifice Spacing = ’,f6.2,’ in’/ 

5 5x,’Orifice Loss Coefficient = ’46.2/) 
pu =s 144. *pu 

psmi = 144.*psmi 
tol = 144*tol 
zliq = zliq/12. 
ztank = ztank/12. 
dsm = dsm/12. 
dsi = dsi/12. 
zsm = zsm/12. 
zsi = zsi/12. 
dzsi = dzsi/12. 
asm = 3.14159*dsm*dsm/4. 
asi = 3.14159*dsi*dsi/4. 
do j = l,nsec 
nl =n2+ 1 
n2 = node(j) 
do i = nl,n2 
as(i) = asec(j)/144. 
enddo 
enddo 

kesm = .5*(1. - (dsi/dsm)**2) 
call value(nb 4 x)d,lod 4 overd 4 overd) 



zu = ztank - zliq 
gc = 32.2 
do 10 j = l,nlim 

qsm = (mdot/asm)**2/(2.*rho*gc) 
call Mct(dsm,ed,mdot,visc,fsm) 
kfsm = fsm* 2 sm/dsm 
kbsm = fsm*loverd 
dpfsm = qsm*(kfsm + kbsm + kcsm) 
dpf 1 = qsm*(kfsm) 
phsm = rho*acc*zsm 
psmo = psmi - dpfsm - phsm 
pback = psmi - dpf 1 
pback = pback/144.0 
dpsm = dpfsm + phsm 
mdoti = mdot/nbar 
qi = (mdoti/asi)**2/(2.*rho*gc) 
call frict(dsi,ed,mdoti,visc/si) 
kbsi = fsi*loverd 
dpsi = qi*kbsi 
pi * psmo - dpsi 
do i = l,n 

x(i) = i*dzsi - dzsi/2. 
aorf = as(i) 
z = ztank - zsi + x(i) 
if (z .It zu) pt = pu 

if (z .gc. zu) pt = pu + rho*acc*(z - zu) 
callpres(i,n,pi,po,pn,ptdpf,ph»mdoti,mdoto,mdots, 
1 dzstdsi,ed,asi,aoif4cs4io,visc,acc) 

if (aorf .le. 0.0) vel(i) « 0.0 
if (aorf .gt 0.0) vel(i) = mdots/(rho*cds*aorf) 
pin(i)=pi 
pOUt(i) =: po 
pnode(i) = pn 
ptank(i) = pt 
mdin(i) » mdod 
mdout(i) = mdoto 
mds(i) - mdots 
pi = pout(i) 
mdod = mdoutd) 
if (i .It n) dpsi = dpsi + dpf - ph 
if (i .eq. n) dpsi * dpsi + dpf - ph/2. 
enddo 

ptcal = pn - ks/(2.*rho*gc)*(mdots/aorf)**2 
delpt(j) = ptcal - pt 
if (j .eq. 1) then 
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if (delpt(j) .It. 0.) mdot = mdot - dmdot 
if (delpt(j) .gt. 0.) mdot = mdot + dmdot 

else 

prod = delpt(j)*delpt(j-l) 

if (delpt(j) .It 0. .and. prod .gt. 0.) mdot = mdot - dmdot 
if (delpt(j) .gt. 0. .and. prod .gt. 0.) mdot = mdot + dmdot 
if (delpt(j) .It 0. .and. prod .It 0.) then 
dmdot = dmdot/2. 
mdot = mdot - dmdot 
endif 

if (delpt(j) .gt 0. .and. prod .It 0.) then 
dmdot = dmdot/2. 
mdot = mdot dmdot 
endif 
endif 

10 continue 

if (abs(delpt(j)) .gt tol) write (6,1002) delpt(j) 

1002 format (’♦♦♦ tank pressure docs not converge, delpt = ', 

1 f8.4,’psi***') 

15 dppump “ (psmi - pt)/144. 
dpsm » dpsm/144. 
dpsi = dpsi/144. 

write (6,200) mdot,dppump,dpsm,dpsi 
200 format (5x,*Pump Flow ^te = ',f6.3,' Ibm/sec'/ 

1 5x,’Pump Pressure Rise = ’,f6.3,' psiV 

2 5x,'Spray Manifold Tube Delta p = *,f6.3,’ psi’/ 

3 5x,'Spray Injection Tube Delta p * *,f6.3,’ psi'//) 
write (6,4) 

4 format (5x,'Node',2x,’Distance’,3x,*Inlet',2x,'Outlef , 

1 3x,'Nodal'3x,Tank’,5x,’Inlef,4x,’Outief, 

2 3x,'Injection',3x,’Injection’,3x,'Orifice') 
write (6,5) 

5 format (24x,'p’,7x,’p',7x,'p’,7x,'p', 

1 7x,'mdot'4x,’mdof,7x,'mdot',6x,'Velocity', 

2 5x,’CdA’) 
write (6,6) 

6 format (13x,'(in)',4x,’(psia)’,2x,'(psia)',2x,'(psia)',2x,'(psia)', 

1 2x,’Gbm/sec)',lx,'Gbm/scc)',lx,’(lbm/sec)',5x,’(^s)', 

2 6x,'(in2)'/) 
do i = l,n 

pin(i) « pin(i)/144. 
poutG) = pout(i)/144. 
pnode(i) = pnode(i)/144. 
ptank(i) * ptank(i)/144. 
as(i) = 144.*as(i) 



x(i) = 12.*x(i) 

write (6,210) i,x(i),pin(i),pout(i),pnode(i),ptank(i), 

1 mdin(i),mdout(i),mds(i),vel(i),as(i) 

210 format (5x,i3,4x,f6.2,3x,4(f6.3,2x), 

1 3(e9.3,lx),4x,f5.2,4x,f7.5) 

enddo 

if (opn .eq. 0.) go to 20 

call crtplt(-l 1,1 1,2, 00,0, l,subhd,xtitl,yp,0.,0., 

1 0.,0.,n,x,pnode,0,0) 

20 if (opt ,eq. 0.) go to 21 

caU crtplt(0 ,0,0,0 1,0,0,0,0,0,0.,0., 

1 0.,0.,n,x,ptank,0,0) 

21 if (omdin .eq. 0.) go to 22 

call crtplt(-l 1,1 1,2,00,0, l,subhd,xtitl,ymdot,0.,0., 

1 0.,0.,n,x,mdin,0,0) 

22 if (omds .eq. 0.) go to 23 

call crtplt(-l 1,1 l,2,00,0,l,subhd,xtitl,ymds,0.,0., 

1 0.,0.,n,x,mds,0,0) 

23 if (ovel .eq. 0.) go to 24 

call crtplt(-l 1,1 l,2,00,0,l,subhd,xtitl,yvel,0.,0., 

1 0.,0.,n,x,vel,0,0) 

24 if (oas .eq. 0.) go to 25 

call crtplt(-l 1,1 1,2,00,0,1, subhd,xtitl,yas,0.,0., 

1 0.,0.,n,x,as,0,0) 

25 return 
end 

subroutine TPHS(HI,RHOL8,RHOG8>lDOT,AN8,HL8JIG8,SL8,SG8, 
1 y8,h8,s8,rho8,p8d) 

double precision hi 4 hol 8 ,rhog 84 il 84 ig 8 ,sl 8 ,sg 8 ,h 8 ,riio 8 ,p 8 d 
real k,mdot 
c 
c 

K=4).41405*(MDOT*MDOT)/(AN8*AN8) 

A1=K*(1.0/RHOG8-1.Q/RHOL8)*(1.0/RHOG8-1.0/RHOL8) 

B1=HG8-HL8+2.0*K/RHOL8*(1.0/RH(X}8-1.0/RHOL8) 

Cl=HL8-HI+K/(RHOL8*RHOL8) 

Y8=(-B1+SQRT(ABS(B1*B1-4.0*A1*C1)))/(2.0*A1) 

IF(Y8.LE.0.00) then 
y8=0.00 
endif 

IF(Y8.GE.1.0) then 
y8=1.00 
endif 

H8=Y8*HG8+(1.0-Y8)*HL8 



S8=Y8*SG8+(1.0-Y8)*SL8 

RHO8=1.0/(1.0/RHOL8+Y8*(L0/RH(Xj8-L0/RHOL8)) 

P8D=144.0/(2.0*32.2*RHO8)*(MDOT/AN8)*(MDOT/AN8) 

return 

END 

subroutine TPS(A,G,HI,HO,K,P,YO,RLO,RGO,YI,RLI,RGI,MDT,PHISQ,b, 
1 d,e,f) 

double precision p 4 :lo 4 ‘go^li/gi»phisq,b,d,e^ 
double precision dpfl,dpf»vo,vi,dpni,dpt,ph,pinlet 
real k,mdt 
c 

RHOO=1.0/(1.0/RLO+YO*(1.0/RGO-1.0/RLO)) 

RHOI=1.0/(1.0/RU+YI*(1.0/RGI-1.0/RLD) 

RHO=(RHOO+RHOD/2.0 
YA=(YI+YO)/2.0 
IF(YA.GT.1.0) then 
ya=0.999 
endif 

RHOLA«(RLI+RLO)/2.0 
aa = ((L0-ya)**2.0)*PfflSQ/(RHOLA) 
bb = LO / rho 
cc = aa 

if (bb .ge. aa) cc = bb 
DPF1«(MDT/A)**2 
DPF=144.0*K*(DPFl)*cc/(2.0*32.2) 

B=DPF 

VO=MDT/(RH(X)*A)*144.0 

VI=MDT/(RHOI*A)*144.0 

DPM=ABS(MDT/A*(VO-VI)/32.2) 

DPT=DPM+DPF 

PH=G*RHO*(HI-HO)/1728.0 

PINLET=P+DPT-PH 

D=PINLET 

E=DPT 

F=DPM 

return 

END 

subroutine value(np^,ypdn,yout) 
c 

c This subroutine performs lagrangian interpolation within a set 
c of (x,y) pairs to give yout coiresponding to xin 
c 

double precision x(np),y(np),xin,yout 


(i-a . 
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if (xin .le. x(l)) yout = y(l) 
if (xin .le. x(l)) return 
if (xin .ge. x(np)) yout = y(np) 
if (xin .gc. x(np)) return 
do 10 i = l,np 
k = i 

if (xin .It x(i)) go to 30 
10 continue 

30 ffr = (xin - x(k - l))/(x(k) - x(k - 1)) 
yout = y(k - 1) + ffr*(y(k) - y(k -1)) 
c write (6,*) np,x,y,xin,yout 
return 
end 

SUBROUTINE VALUE2(NPX,NPY,X,Y,Z,XIN,YIN^OUT) 
c DIMENSION X(NPX),Y(NPY),Z(NPX,NPY) 

double precision X(NPX),Y(NPY)^(NPX,NPY)^IN,YIN^OUT 
C 

DFCXIN .LE. X(l) .AND. YIN .LE. Y(l))ZOUT=Z(l,l) 

JFQW .LE. X(l) .AND. YIN .LE. Y(1))RETURN 
C 

IF(XIN .GE. X(NPX) .AND. YIN .LE. Y(l))ZOUT=Z(NPX,l) 

IF(XIN .GE. X(NPX) .AND. YIN .LE. Y(1))RETURN 
C 

IF(XIN .LE. X(l) .AND. YIN .GE. Y(NPY))ZOUT=Z(l,NPY) 

IF(XIN .LE. X(l) .AND. YIN .GE. Y(NPY))RETURN 
C 

W(xm .GE. X(NPX) .AND. YIN .GE. Y(NPY))ZOUT=Z(NPX,NPY) 
worn .GE. X(NPX) .AND. YIN .GE. Y(NPY))RETURN 
C 

wpw .GT. X(l))GO TO 30 
DO 20 1=1 ^PY 
M=I 

IF(YIN .LT. Y(I))GO TO 25 
20 CONTINUE 

25 FFRY=(YIN-Y(M-1))/(Y(M)-Y(M-1)) 
ZOUT=Z(l,M-l)+FFRY*(Z(lM)-Z(lM-l)) 

RETURN 

C 

30 IF(XIN .LT. X(NPX))GO TO 60 
DO 50I=1J^PY 
M=I 

IF(YIN .LT. Y(I))GO TO 55 
50 CONTINUE 
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55 FFRY=(YIN-Y(M-1))/(Y(M)-Y(M-1)) 

ZOUT=Z(NPX,M- 1)+FFRY*(Z(NPX,M)-Z(NPX,M- 1)) 

RETURN 

C 

60 IF(YIN .GT. Y(l))GO TO 90 
DO 80 1=1, NPX 
L=I 

IF(XIN XT. X(I))GO TO 85 
80 CONTINUE 

85 FFRX=(XIN-X(L-1))/(X(L)-X(L-1)) 
ZOUT=Z(L-l,l)+FFRX*(Z(L,l)-Z(L-l,l)) 

RETURN 

C 

90 IF(YIN .LT. Y(NPY))GO TO 120 
DO110I=l,NPX 
L=I 

IF(XIN .LT. X(I))GO TO 1 15 
no CONTINUE 

115 FFRX=(XIN-X(L-1))/(X(L)-X(L-1)) 
ZOUT=Z(L-l,NPY)+FFRX*(Z(L,NPY)-Z(L-l,NPY)) 

RETURN 

C 

120 DO 130I=1,NPX 
L=I 

W(xm .LT. X(I))GO TO 135 
130 CONTINUE 

135 FXR=(XIN-X(L-1))/(X(L)-X(L-1)) 

DO 140 1=1 J^PY 
M=I 

IF(YIN .LT. Y(I))GO TO 145 
140 CONTINUE 

145 FYR=(YIN-Y(M-l))/Or(M)-Y(M-l)) 

C 

ZXLO=Z(Ln,M-l)+FYR*(Z(L-lJ4)-Z(L-l>I-l)) 

ZXHI=Z(L>!-1)+FYR*(Z(L,M)-Z(L>!-1)) 

C 

ZOUT=ZXLO+FXR*(ZXHI-ZXLO) 

RETURN 

END 

SUBROUTINE VALUE3(NPX,MAXYJ^PY,X,Y,Z;aN,YINZOUD 
integer npy(maxy) 

double precision X(NPX),Y(NPX>IAXY),Z(NPX,maxy),xin,yin,zout 
c double precision X(15),Y(15,15),Z(15,15),xin,yin,zout 
dimension a(50),b(50),c(50),d(50) 



c 


IF(XIN .LE. X(l) .AND. YIN .LE. Y(l,l)) then 
Z0UT=Z(1,1) 

RETURN 
, endif 
C 

IF(XIN .GE. X(NPX) .AND. YIN .LE. Y(npx,l)) then 
ZOUT=Z(NPX,l) 

RETURN 

endif 

C 

IFQW .LE. X(l) .AND. YIN .GE. Y(1,NPY(1))) then 
Z0UT=Z(1^PY(1)) 

RETURN 

endif 

C 

IF(XIN .GE. X(NPX) .AND. YIN .GE. Y(npx,NPY(npx))) then 
ZOUT=Z(NPX^Y(npx)) 

RETURN 

endif 

C 

IF(XIN .GT. X(l))GO TO 30 
DO I*1,NPY(1) 

M=I 

IF(YIN .LT. Y(14))GO TO 25 
enddo 

25 FFRY=(YIN-Y(1,M-1))/(Y(1,M)-Y(1,M-1)) 
ZOUT^Z(l,M-l)+FFRY*(Z(l>l>Z(l,M-l)) 

RETURN 

C 

30 JFQW .LT. X(NPX))GO TO 60 
DOI=l.NPY(npx) 

M=I 

IF(YIN .LT. Y(npx,I))GO TO 55 
enddo 

55 FFRY=(YIN-Y(npx,M-l))/(Y(npx,M)-Y(npx,M-l)) 
ZOUT=Z(NPX,M-l)+FFRY*(Z(NPX>l)-Z(NPX>I-l)) 
RETURN 
C 

60 DOI=lJIPX 
L=I 

IF(XIN .LT. X(I))GO TO 85 
enddo 

85 doj = l,npy(l) 

a(i)=yGd) 
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b(j)=zaj) 

enddo 

call value(npy(l),a,b,yin,zl) 
do k = 1, npy(l-l) 
c(k) = yG-14c) 
d(k) = z(l-l,k) 
enddo 

call value(npy(l-l),c,d»yin,z2) 

FFRX=(XIN-X(L- 1))/(X(L)-X(L- 1)) 

Z0UT=Z2+FFRX*(Z1-Z2) 

RETURN 

END 

subroutine XPARAM(YO,VLO,VGO.RLO,RGO,YI,VLI,VGI,RLI,RGI,x) 
c 

double precision vlo,vg 04 ’lo,rgo,vli,vgi,rli,rgi,x 
double precision visla,visga,rhola^hoga,xl 
c 

YA:=(YI+YO)/2.0 
- VISLA=(VU+VLX))/2.0 
VISGA*(VGI+VGO)/2.0 
RHOGA=(RGI+RGO)/2.0 
RHOLA=(RLI+RLO)/2.0 
IF(YA.LE.1.0E-06)then 
ya*1.0e-06 
end^ 

IF(YA.GE.1.0) then 
ya=1.0 
endif 

X1=((ABS(VISLAATSGA))**0.1)*(((1.0-YA)/YA)**0.9) 

X= Xl*SQRT(ABS(RHOGA/RHOLA)) 
c write (6,10) yi,yo,ya,vli,vlo,visla,vgi,vgo,visga 4 ‘gi,rgo 4 hoga, 
c 1 rli 4 'lo,rhoia,xl,x 

c 10 format (6(3x,f8.4)y,6(3x.f8.4)y.3(3x,f8.4),2(3x,f 10.4)) 
return 
END 

3.4.2 Integrated Z ero-g TVS Model 
c 

c program tvs.f to model zero-g Thermodynamic Venting System (TVS) 
c transient performance 
c 

character* 1 label 

real ptitl( 1 8,4),subhd( 1 8),xtitl( 1 8)4cetitl(4,2), 

1 ypu(18),ytu(18),ymu(18),ypl(18),ytl(18),yml(18), 



2 ytw( 1 8) ,y mwl( 1 8),ymdv( 1 8),ymdsu( 1 8),ymddu( 1 8),ymdb w( 18), 

3 ymdif(l 8),ymds(l 8),ynpump(l 8),ydppmp( 1 8) 
real tsat(25),psat(25),enthf(25),enthg(25), 

1 shpf(25),densf(25),condf(25),viscf(25),texpf(25) 
real pvap(8),tvap(8),tnrm(llXenth(8,ll),shv(8,ll), 

1 shp(8, 1 l),dens(8, 1 l),cond(8,l l),visc(8, 1 1) 

real tal( 1 2),s hpal( 1 2),node( 1 0),dorf ( 1 0) 

real as ( 1 00),ddrop( 1 00) ,ad( 1 00), vd( 100) ,veld( 1 00),ndrop( 1 00) 

real ton( 1 00), toff ( 1 00),tcyc( 1 00),pton( 100), 

1 mvent(100),mvnt(100),sfc(100) 

real mup,mlp,mwlp,mvp 4 ndvp 

real mdsp,mdslp,mdsup,mddup,indbwp,mdlup,npumpp 
real mu 4 iil,mwl,mw,mv,miu,rnii 4 nwli,mvi 4 nwlrTiax 
real kwu,muwu4cuwl,muuwl,kwl,muwl,muul4cul 
real kus,muus Jds,muls Jmd,muud,nud 

real ks,mds,mdsl,mdsu,iiiddu,mdsw,mdbw,mdbwnix,mdlu,mdul 
realmdv,mdvent,mdp,mdpuinp,mdsne,mdsul,mdcond,mdcndp,mdsi 
real rndotd,npumpd,npump,npumpi,nmax,nrnin,ip 
integer jsymb(2) 
integer*2 ibox41oc 

conimon/pltcoin/niisc(3),nc,miss(13),dclini,ltick,nfig,nptmin, 

1 nUnes,nchlin,ptitl 

common /contrlAboxJloc 

common /sprayin/dsm,zsm,dsi,zsi,norf,nbar, 

1 ks,cds,as,vd,roverd,dmdot,tol,nlim, 

2 rhos,viscs 

common /pumpin/npumpd,nmin,nmax,hpid,hpmotr,effp,ip,dm 

common /tankdim/dtank,hcyl,hbulk,tkw,cl 

common /tankvoWl,v2,vbi,vbo,vci,vco 

common /tankarea/ab,ac,at 

common /tankout/j4iline4print,tiinep(3000), 

1 pup(3000),tup(3000),vup(3000),mup(3000),twp(3000), 

2 plp(3000),tlp(3000),vlp(3000),mlp(3000), 

3 pwlp(3000),twlp(3000),vwlp(3000),mwlp(3000), 

4 mvp(3000),tsp(3000),mdlup(3000),mdsp(3000),mdslp(3000), 

5 mdsup(3000),mddup(3000)4ndbwp(3000),mdulp(3000),mdcndp(3000), 

6 qwup(3000),quwlp(3000),qwlp(3000),qulp(3000),qudp(3000), 

7 qusp(3000),qlsp(3000),npumpp(3000),dppmp(3000) 
nc = 2 

nfig = 0 
nlines = 4 
ibox = 257 
jsymb(l) ss 40 
jsymb(2) = 43 

dataketitl/4hEVAP,4hORAT,4hTION,4h , 



1 4hCOND,4hENSA,4hTION,4h / 

data ni,rhow,gc,ed/766.0, 176.3, 32.2, l.Oe-6/ 

data tal/0.0,25.0,50.0,75.0,100.0,125.0, 150.0,180.0,260.0, 

1 485.0,760.0, lOOO.O/nal/12/ 

data shpal/0.0,0.001,0.0075, 0.02,0.0525,0.09,0.11,0.125,0.15, 
1 0.2,0.225,0.25/ 


c 

c read input data for tank model 
c 


read (5,100) label 

read ( 5 *) xd,xchar,he,xcond,prtsp,outp 
read (5,100) label 

read (5,*) pui,tui,pli,twi,twli,fuU,xll 
read (5,100) label 

read (5,*) dtank,hcyl,hbulk,tkw,dsb,dl,d2 
read (5,100) label 

read (5,*) mdvent,mdsi,dthex,qflux,g,hwliq 
read (5,100) label 

read (5,*) pmin,pmax,delt2,ipmt2,iplot2 
read (5,100) label 

read (5,*) fintim,deltl,xdeltl4pmtl4plotl,nline 

read (5,100) label 

read (5,*) opu,otu,omu,opl,otl,oml 

read (5,100) label 

read (5,*) otw,omwl,omdv,omdsu,omddu,omdbw 
read (5,100) label 

read (5,*) onxUu,onidul,onxls,onpump,odppmp 

read (5,200) subhdpctitl,ypu,ytu,ymu,ypl 

read (5,200) ytl,yml,ytw,ymwl,ymdv,ymdsu 

read (54^00) ymddu,ymdbw,ymdif,ymds,ynpump,ydppmp 

read (5,201) ptitl 


c 

c read input data for pump model 
c 


read (5,100) label 

read (5,*) mdotd,dp(M}pumpi,npumpd,xhp,xn 
read (5,100) label 
read (5,*) deltat,e^ 


c 

c read input data for spray manifold^jection tubes model 
c 


read (5,100) label 

read (5,*) dsm,zsm,dsi,zsi 4 iorf,nbar 
read (5,100) label 

read (5,*) ks,cds,roverd,dmdot,tol,nlim 






m-58 



read (5,*) nsec 
do i = l,nsec 

read (5,*) node(i),dorf(i) 
enddo 

100 format (//al) 

200 format (18a4/18a4/18a4/18a4/18a4/18a4) 

201 format (18a4) 
c 

c read LH2 saturation properties 
c 

open(unit=2,fUe=’h2prop’,status='old') 
read (2,*) nsat 
do i = l,nsat 

read (2,*) tsat(i),psat(i),enthf(i),dunivar,shpf(i),densf(i), 

1 texpf(i),condf(i),viscf(i) 

read (2,*) enthg(i) 
enddo 
c 

c read GH2 properties as a function of pressure and temperature 
c 

read (2,*) np,nt 

read (2,*) (tnrm(i),i=l,nt),tcnst 

do i = l,np 

read (2,*) pvap(i),tvap(i) 
do j = l,nt 

read (2,*) enth(ij),shv(i j),shp(io),dens(i j), 

1 cond(i,j),visc(ij) 

enddo 
enddo 

if (outp .eq. 1.0) then 
write (6,1) 

1 format (5x,TANK DIMENSIONS’/) 
write (6,2) dtank,tkw,hcyl,hbulk 

2 format (5x,Tank Diameter ',f6.1,'in'/ 

1 5x,Tank Wall Thickness ',f6.2,’ in'/ 

2 5x,’Cylinder Height ’,f6.1,' in’/ 

3 5x,’Bulkhead Height ',f6.1,' in'//) 

write (6,3) 

3 format (5x,’SPRAYMANIFOLD/INJECnON TUBE DIMENSIONS’/) 
write (6,4) zsm,dsm,zsi,dsi,nbar,norf,ks,cds 

4 format (5x,’Spray Manifold Tube Length ',f6. 1 ,' in'/ 

1 5x,’Spray Manifold Tube ID '4^6.3,' in'/ 

2 5x,’Spray Injection Tube Length '46. 1 ,' in’/ 

3 5x,'Spray Injection Tube ID ',f6.3,' in'/ 

4 5x,’Number of Spray Injection Tubes '46/ 
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5 5x, 'Number of Orifices ',i6/ 

6 5x, 'Orifice Loss Coefficient ',f6.2/ 

7 5x, 'Orifice Discharge Coefficient ',f6.2/) | 

write (6,5) 

5 format (t2,'Time',tl0,'pU',tl8,'TU’,t26,’VU',t34,'MU’, 

1 t42,'pL',t50,’TL',t58,’VU,t66,'ML',t74,'TW, | ; 

2 t82,'TWL’,t90,'Npump’,t98,'dppump',tl06,'dTpump’, ‘ 

3 tl 14,'HPO',tl22,'mdS',tl30,'mdSU’,tl38,'mdDU', 

4 tl46,'mdBW,tl54,'mdLU',tl62,'mdcond') | 

write (6,6) 

6 format (t2,'sec',t9,'psia',tl9,'R’,t25,'ft3’,t33,’lbm', r 

1 t41,'psia',t51,’R',t57,’ft3',t65,’lbm’,t75,'R’, 

2 t83,’R’,t9i;rpm',t99,'psid',tl09.’R', 

3 tl 14,'HF,tl20,'lbm/sec',tl28,'lbm/sec’,tl36,’lbm/sec', 

4 tl44,'lbm/sec’,tl52,'lbm/sec’,tl60,'lbm/secy) 
endif 

do j = l,nsec f 

nodel = node2 + 1 I 

node2 ^ node(j) 
do i - nodel 4 iode 2 
ds = dorf(j)/12.0 
as(i)=3.14159*ds**2/4.0 
ddrop(i) = xd*ds 
ad(i) = 3.14159*ddropa)**2 
vd(i) = 3.14159*ddrop(i)**3/6.0 

enddo | 

enddo ^ 

dtank - dtank/12.0 
hcyl = hcyl/12.0 
hbulk = hbulk/12.0 

htank = hcyl + 2,0*hbulk ^ 

tkw = tkw/12.0 

dsb = dsb/12.0 

dsm = dsm/12.0 

zsm = zsm/12.0 

dsi » dsi/12.0 

zsi = zsi/12.0 

xll=xll/12.0 

dl=dl/12.0 

d2 = d2/12.0 

dzsi = zsi/norf | 

debar - xchar^dtank 

C 5.' 

tv 

c initial ullage and liquid masses | 

c 

I 
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call volarea(vt) 

dtanki = dtank - 2.0*tkw 

call value(nsat,psat,tsat,pli,tli) 

call value(nsat,tsat,densf,tli^hol) 

call value(nsat,tsat,densf,twli 4 ‘howl) 

vli = full/1 00.0*vt 

mli = rhol*vli 

vui = vt - vli 

mui = 144.0*pui*vui/(ni*tui) 
ts = tli 
tsw = ts 


c 

c in 1 g, set maximum thickness of wall liquid layer to 0.01 in 
c due to liquid run-off and calculate maximum wall liquid mass 
c 

call area(vli,y»awu,aul,awl,hliq,hu) 
if (g .ge. 1.0) mwlmax = rhol*awu*0,01/12.0 
c 

c pump design conditions 
c 

dm = 720.0/(3. 14159*npumpd)*(2.0*144.0*dpd/rhol)**0.5 

hpid = xhp*mdotd*dpd*144.0/(550.0*ef^*rhol) 

hpmotr = hpid 

nmax = npumpd*(1.0 + xn) 

nmin = npumpd*(1.0 - xn) 

ip = 6.018e+05*hpid*deltat/(npumpd*(nmax - nmin))/2.0 
if (pui .ge. pmax) mdp » mdsi 
if (pui .le. pmin) mdp = 0.0 
mdpump = mdsi 
c 

c time integration of variables 
c 

i = 0 
j = 0 

85 if (npump .gt 0.0) then 
delt deltl 

if (pu .le. pi) delt = xdeltl^deltl 
iprint = ipmtl/xdeltl 
iplot » iplotl 
else 

delt = delt2 
iprint =* ipmt2 
iplot = iplot2 
endif 

mu = mui + dmudt*delt 
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ml = mli + dmldt*delt 
if (ml .le. 0.0) mi = 0.0 
mwl = mwli + dmwldt*delt 
if (mwl .le. 0.0) mwl = 0.0 

if (g .ge. 1.0 .and. mwl .ge. mwlmax) mwl = mwlmax 

mv = mvi + mdv*delt 

vl = vli + dvldt*delt 

if(vl .le. 0.0) vl = 0.0 

vwl = mwl/rhowl 

vu = vt - vl - vwl 

tl = tli + dtidt*delt 

tu = tui + dtudt*delt 

call value(nsat,psat,tsat,pu,tusat) 

if (tu .le. tusat) tu = tusat 

twl = twli + dtwldt*delt 

if (twl .ge. tusat) twl = tusat 

if (twl .le. ts) twl =2 ts 

tw = twi + dtwdt*delt 

qpump = qpumpi + hpo*0.707*delt 

mui = mu 

mli = ml 

mwli = mwl 

mvi = mv 

vli = vl 

tui = tu 

tli = tl 

twli = twl 

twi = tw 

qpumpi = qpump 


c ullage, bulk liquid, and wall liquid pressures 
c 

pu = mu*ru*tu/(144.0*vu) 
call value(nsat,tsat,psat,tl,pl) 
call value(nsat,tsat,psat,twl,pwl) 
c 

c pump control logic 
c 

if (pu .ge. pmax) flagl = 1.0 
if (pu .le. pmin) flagl = 0.0 
c 

c performance calculations 
c 

if (flagl .eq. 0.0) then 
if (eye .cq. 0.0) ncyc = ncyc + 1 
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toff(ncyc) = time - timeon 
timeof = time 
eye = 1.0 
endif 

if (flagl .eq. 1.0) then 

ton(ncyc) = time - timeof 
tcyc(ncyc) = ton(ncyc) + toff(ncyc) 
pton(ncyc) = 100.0*ton(ncyc)/tcyc(ncyc) 
sfc(ncyc) = 3600.0*mv/time 
mvent(ncyc) = mv 
pup(ncyc) = pu 
tup(ncyc) = tu 
mup(ncyc) = mu 
plp(ncyc) = pi 
tlp(ncyc) = tl 
mlp(ncyc) = ml 
twp(ncyc) = tw 
mdsp(ncyc) = mds 
timeon = time 
eye = 0.0 
endif 
e 

e pump model 
e 

eall value(nsat,tsat,shpf,tl,cpl) 
nmdot = 0 

87 eall pump(flagl 4 ipump^pumpi 4 elt,mdpump^hol,epl,dpptimp,dtpump, 
1 dnt,hpo) 

if (dppump .le. 0.01) mdp = 0.0 
if (dppump .gt 0.01) mdp = mdpump 
tpump = tl + dtpump 
e 

e vent eontrol logic 
c 

psatp = -36.37 + 4.6054*tpump - 0.20369*tpump*tpump 
1 + 0.003 1745*tpump*tpump*tpump 

if (psatp .gt pmin) flag2 = 1.0 
if (psatp .Ic. pmin) flag2 = 0.0 
mdv = 0.0 

if (flagl .gt 0.0 .and. flag2 .gt 0.0) mdv = mdvent 
c 

c pressure drop in the ^circulation line (between the pump outlet 
c and spray manifold inlet) 


if (mdp .Ic. 0.0) dprec = 0.0 



if (mdp ,gt 0.0 ) then 
call value(nsat,tsat4ensf,tpump»rhop) 
call value(nsat,tsat,viscf,tpump,viscp) 
call pdrop(mdp4'hop,viscp,xll,dl»d2,dprec) 
endif 
c 

c spray manifoldAnjection tube model 
c 

call aiea(vl,y,awu,aul,awl»hliq,hu) 
dhn htank - hu - dzsi/2.0 
ptn * pu + rhol*g*dhn/144,0 
psmi « ptn + dppump - dpiec 
if (mdp .gt. 0.0) ts = tpump • dthex 
if (mdp .le. 0.0) ts = tpump 
call value(nsat»tsat,densf,ts^hos) 
call value(nsat,tsat,viscf,ts, vises) 
if (mdp .gt 0.0) then 

c2jlspray(psmi,mdp,pu»hu4iliq,htanltdchar/hol,g,ed, 
1 time,prtsp,mds4Tidsu4ndsl,veld,ndrop,norfu) 
if (mds .le. 0.0) mds » mdsi 
mdpump * (mdp + mds)/2.0 
delmd = mdp - mds 
if (abs(delmd) .It 0.01) go to 86 
nmdot = nmdot + 1 
if (nmdot .It 10) go to 87 
write (6,8) time,delmd 

8 format (’*♦♦ pump flow rate does not converge at time 
1 flO.2/ sec, delmd « '48.4,' Ibm/sec ***') 
else 

mdsu » 0.0 
mdsl » 0.0 
endif 
c 

c heat*transfer rates 
c 

c wall-to-uUage 
c 

86 if (mwl .gt 0.0) qwu s 0.0 
if (mwl .le. 0.0) then 
twu = (tw + tu)/2.0 
betawu = 1.0/twu 
twun = (twu - tusat)/(tcnst - tusat) 
call value2(np,ntpvap,tnrm,shp,pu,twun,cpwu) 
call valuc2(np,nt,pvap,<nrm,dens,pu,twunihowu) 
call value2(np,nt,pvap,tnrm,cond,pu,twun4cwu) 


call value2(np,nt,pvap,tnrm,visc>pu,twiin,muwu) 
call htc(tw,tu,dtanki,cpwu 4 'howu,betawu,kwu,muwu,g,hwu) 
qwu = hwu*awu*(tw - tu) 
endif 

iillage-to-wall liquid 

if (mwl .le. 0.0) quwl = 0.0 
if (mwl .gt 0.0) then 
tuwl =s (tu + twl)/2.0 
betuwl = 1.0/tuwl 
call value(nsat,psat,tsat,pu,tusat) 
tuwln “ (tuwl - tusat)/(tcnst - tusat) 
call value2(np»nt,pvap,tnrm,shp,pu,tuwln,cpuwl) 
call value 2 (np,nt,pvap,tnrm,dens,pu,tuwln 4 ‘houwl) 
call value2(np,nt,pvap,tnrm,cond,pu»tuwln4cuwl) 
call value2(np»nt,pvap,tnim,visc,pu,tuwln,muuwl) 
call htc(tu,twl 4 tanki,cpuwl 4 iiouwl,betuwl 4 aiwl,muuwl»g,huwl) 
quwl = huwl*awu*(tu - twl) 
endif 

wall-to-wall liquid 


if (mwl .le. 0.0) qwl = 0.0 
if (mwl .gt 0.0) then 
call value(nsattsat,shpf,twl,cpwl) 
call value(nsattsat,texpf,twl,betawl) 
call value(nsat,tsat,condf,twl4cwl) 
call value(nsattsatviscf,twl 4 nuwl) 
call value(nsat,tsat,densf,twUhowl) 
hwl = 0.13*(3600.0*32.2*g*betawl*abs(tw - twl)*ihowl**2 
1 *kwl**2*cpwl/muwl)**(1.0/3.0) 

if (hwl .gt hwliq) hwl = hwliq 
qwl = hwl*awu*(tw - twl) 
endif 

uUage-to-liquid 

tul = (tu + tl)/2.0 
betaul = 1.0/tul 

tuln = (tul - tusat)/(tcnst - tusat) 
call value 2 (np 4 itpvap,tnnn,shp»pu,tuln,cpul) 
call value2(np»ntpvap,tnim,dens,pu,tuln^houl) 
call value2(np,nt,pvap,tnrm,cond,pu,tuln»kul) 
call value2(np^tpvap,tnrm,visc,pu,tuln,muul) 



callhtc(tu,tl4tanki,cpul,rhoul,betaul,kul,muul,g,hul) 
qul = hul*aul*(tu - tl) 
c 

c ullage-to-droplet 
c 

td = ts 

tud = (tu + td)/2.0 

tudn ss (tud - tusat)/(tcnst - tusat) 

call value2(np,nt,pvap,tnrm,dens,pu,tudn4ioud) 

call value2(np,nt,pvap,tnrm,visc,pu,tudn,muud) 

call value2(np,nt,pvap,tnnn,cond,pu,tudn4cud) 

qud = 0.0 

if (mdp .le. 0.0) go to 17 
do io = l,norfu 

red = rhoud*veld(io)*ddrop(io)/muud 
c 

c textbook heat>transfer correlation for liquid droplets (Kreith) 
c 

nud = 0.3125*red**0.602 
if (ddrop(io) .gt 0.0) hud - nud*kud/ddrop(io) 
qud = qud + ndrop(io)*hud*ad(io)*(tu - td) 
enddo 

17 call value(nsat,tsat4ensf,td^od) 


c 

c ullage-to-spray bar 
c 

tus = (tu + ts)/2.0 
betaus = 1.0/tus 


tusn ~ (tus • tusat)/(tcnst - tusat) 
call value 2 (np 4 it,pvap,tnrm,shp,pu,tusn,cpus) 
call value2(np,nt,pvap,tnrm,dens,pu,tusn,rhous) 
call value 2 (np 4 it,pvap,tnnn,cond,pu,tusn 4 ois) 
call value 2 (np 4 it,pvap,tnim,visc»pu,tusn,muus) 
call htc(tu,ts,dsb,cpus,rhous,betaus4nis,muus ,g,hus) 
aus = 3.14159*dsb*hu 


qus = nbar*hus*aus*(tu - ts) 
qus = 0.0 


c 


c liquid-to-spray bar 


c 


tls = (tl + ts)/2.0 
call value(nsat,tsat,shpf,tls,cpls) 
call value(nsat,tsat,densf,tls»ihols) 
call value(nsat,tsat,texpf4s,betals) 
call value(nsat,tsat,condf,tls^s) 
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call value(nsat,tsat,viscf,tls, mills) 
call htc(tl,ts,dsb,cplsjhols,betals4ds,muls,g, his) 
als = 3.14159*dsb*hliq 
qls = nbar*hls*als*(tl - ts) 
qls = 0.0 - 

c 

c environment-to-wall 
c 

qew = qflux*awu 
qel = qflux*awl 
c 

c mass-transferrates 
c 

c droplet-to-ullage boil-off 
c 

call value(nsat,psat,enthf,pu,hf) 
call value(nsat,psat,enthg,pu,hgsat) 
hfgu = hgsat - hf 

mddu = (qud/3600.0 - mdsu*cpl*(tusat - ts))/hfgu 
if (mddu .It 0.0) mddu = 0.0 
if (mddu .It 0.0) qud = 3600.0*mdsu*cpl*(tusat - ts) 
if (mddu .It 0.0) ts = tusat - qud/(3600.0*mdsu*cpl) 
if (mddu .gt mdsu) mddu - mdsu 

if (mddu .gt mdsu) qud = mdsu*(hfgu + cpl*(tusat - ts))*3600.0 
c 

c non-evaporated spray droplet 
c 

mdsne = mdsu - mddu 
if (g .ge. 1.0) then 
mdsw = 0.0 
mdsul s mdsne 
else 

mdsw = mdsne 
mdsul 0.0 
endif 
c 

c droplet boil-off from wall 
c 

call value(nsat,tsatenthf,tl,hf) 
call value(nsattsat,enthg,tl4ig) 
hfgl = hg - hf 

dpudt = pu*(dmudt/mu + dtudt/tu - dvudt/vu) 
if (mwl .le. 0.0) mdbw = 0.0 
if (mwl .gt 0.0) then 
if (pu ,gt pwl) mdbw = 0.0 


\ 

I • 


if (pu .It pwl .and. dpudt .gt 0.0) 

1 mdbw = ((qwl + quwl)/3600.0 - mdsw*cpl*(twl - tsw))/hfgl 
if (pu .It pwl .and. dpudt .le. 0.) then 
dtdp = 0.37781 - 4.9170e-3*pwl + 21.7623e-6*pwl*pwl 
mdbw = ((qwl -H quwl)/3600.0 - mdsw*cpl*(twl - tsw) 

1 - mwl*cpl*dtdp*dpudt)/hfgl 

endif 
endif 

mdbwmx = mwl/delt 

if (mdbw .ge. mdbwmx) mdbw = mdbwmx 


c 

c Uquid-to-uUage boil-off 
c 


if (pu .gt pi) mdlu = 0.0 
if (pu .le. pi .and. dpudt .gt 0.0) 

1 mdlu = (ql/3600.0 - mdsl*cpl*(tl - ts))/hfgl 
if (pu .le. pi .and. dpudt .le. 0.0) then 
dtdp = 0.37781 - 4.9170e-3*pl + 21.7623e-6*plV 
mdlu = (ql/3600.0 - mdsl*q>l*(tl - ts) 

1 - ml*cpl*dtdp*dpudt)^fgl 

endif 


c 

c condensation 
c 

if (tu .gt tusat) mdul = 0.0 

if (tu .le. tusat) mdul = (qud + qul + quwl)/(3600.0*hfgu) 
if (he .eq. 0.0 .and. flagl .eq. 1.0) 

1 mdcond = xcond*qul/(3600.0*hfgu) 
if (flagl .eq. 0.0) mdcond = 0.0 
if (he .eq. 1.0) mdcond = 0.0 
c 

c rates of change of ullage, wall liquid, and bulk liquid masses 
c 

dmudt - mddu + nxlbw ■i' mdlu - mdul - mdcond 
dmwldt = mdsw > mdbw 

dmldt * mdsl + mdul + mdcond + mdsul - mdlu - mds - mdv 
if (g .ge. 1.0 .and. mwl .ge. mwlmax) dmldt = dmldt + dmwldt 
c 

c rates of change of ullage, liquid, and wall liquid volumes 
c 

call value(nsattsatdensf,tl,rhol) 

dvldt == dmldt^hol 

if (mwl .le. 0.0) dvwldt = 0.0 

if (mwl .gt 0.0) dvwldt = dmwldt/rhowl 

dvudt = -dvldt - dvwldt 


f. 




r- 

I: 
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rates of change of temperature 
ullage 

qu = qwu - quwl - qul - qud - qus 

if (mdcond .eq. 0.0) qu = qwu - quwl - qud - qus 

enthu = dmudt*hgsat 

tun = (tu - tusat)/(tcnst - tusat) 

call value2(np,nt,pvap,tnrm,shv,pu,tun,cvu) 

dtudt = (qu/3600.0 - 144.0/778.0*pu*dvudt + enthu 

1 - cvu*tu*dmudt)/(mu*cvu) 

bulk liquid 

ql = qel + qul - qls 

if (mddu .It 0.0) td = ts 

if (mddu .ge. 0.0) td = tusat 

if (ml .gt 0.0) dtldt = (ql/3600.0 - mdlu*hfgl 

1 + mdsul*cpl*(td - d) - mdsl*cpl*(tl - ts))/(ml^cpl) 

if (ml .le. 0.0) dddt =* 0.0 

wall liquid 

qudchk = 3600.0*mdsu*cpl*(tusat - ts) 
if (qud .gt qudchk) tsw = tusat 

if (qud .le. qudchk) tsw = ts + qud/3600.0/((mdsu + 0.000 l)*cpl) 
if (mwl .gt 0.0) dtwldt = ((qwl + quwl)/3600.0 
1 - mdsw*cpl*(twl - tsw))/(mwl*cpl) 

if (mwl .le. 0.0) dtwldt = 0.0 

tank wall 

call vwaU(hliq,vw) 
mw = rhow*vw 
call value(nal,tal,shpal«tw,cpw) 
qw = qew - qwu - qwl 
dtwdt = qw/3600.0/(mw*cpw) 

output listing 

if (outp .eq. 0.0) go to 19 
if (mod(i4print) .eq. 0.0) 

1 write (6,7) timc,pu,tu,vu,mu,pl,d,vl,ml,tw,ts,npump, 

2 dppump,dtpump,hpo,mds,mdsu,mddu,mdbw,mdlu,mdcond 



7 format (f7.0,f6.2,9(2x^6.2),2x,f6.1, 
1 3(2x,f6.3),6(lx/7.4)) 

if (mod(i,iplot) .ne. 0.0) go to 19 
j=j + l 
timep(j) = time 
pup® = pu 
tup(j) = tu 
vup(j) = vu 
n^up(j) = mu 
twp(j) = tw 
plp(i) = pl 
tlpO) - tl 
vlp® = vl 
mlp(j) = ml 
pwlp® = pwl 
twlp(j) = twl 
vwlpQ) = vwl 
mwlp(j) = mwl 
mvp® = mv 
i™lvp(j) = mdv 
tsp(j) * ts 
mdlup(j) = mdlu 
mdsp(j) = mds 
mdslp(j) = mdsl 
nulsup(j) = mdsu 
mddup® = mddu 
mdbwp(j) = mdbw 
mdulp(j) = mdul 
mdcndp(j) = mdcond 
qwup(j) = qwu 
quwlp(j) = quwl 
qwlp® qwl 
qulp(j) = qul 
qudp(j) = qud 
qusp(j) = qus 
qlsp(j) = qls 
npumpp(j) = npump 
dppmp(j) = dppump 
19 i = i+l 

time - time + delt 
if (time .le. fintim) go to 85 
if (ton(ncyc) .eq. 0.0) ncyc = ncyc - 1 
write (6,999) 

999 format Cl’) 
do i = l^icyc 



ston = ston + ton(i) 
stoff = stoff + toff(i) 
stcyc = stcyc + tcyc(i) 
if (i .eq. 1) mvnt(i) = mvent(i) 
if (i .gt 1) mvnt(i) = mvent(i) - mvent(i-l) 
write (6,1000) i,ton(i),toff(i),tcyc(i),pton(i),mvnt(i),sfc(i) 

1000 format (5x, 'Cycle No. = ',il2/ 

1 5x,’On Time = ',fl2.3,' sec’/ 

2 5x,’Off Time = ',fl2.3,' sec'/ 

3 5x,’Cycle Time * ',fl2.3,' sec’/ 

4 5x,’%OnTime = ’,fl2.3,’ %'/ 

5 5x,'Vented Mass = ',fl2.3,' Ibm’/ 

6 5x,'Specific Fuel Consumption = ',f 12.3,' Ibnt/hr'/) 
write (6,1001) pup(i),tup(i),mup(i),plp(i),tlp(i),mlp(i), 

1 twp(i),mdsp(i) 

1001 format (5x, 'Ullage pressure = *,fl2.3,' psia'/ 

1 5x,'Ullage temperature = ',fl2.3,' R'/ 

2 5x,'Ullage mass = ',fl2.3,' Ibm’/ 

3 5x,’Liquid pressure = ',fl2.3,' psia'/ 

4 5x,'Liquid temperature = ',fl2.3,' R'/ 

5 5x,'Liquid mass = ',fl2.3,' Ibm'/ 

6 5x,'Wall temperature = ',fl2.3,' R'/ 

7 5x,'Pump flow rate = ',fl2.3,' Ibm/sec’/) 

enddo 

ptonav = 100.0*ston/stcyc 

sfcav = 3600.0*mvent(ncyc)/stcyc 

write (6,1002) ston,stoff,stcyc,ptonav,mvent(ncyc),sfcav 

1002 format (5x,'On Time = ',fl2.3,' sec'/ 

1 5x,'Off Time = ',fl2.3,' sec’/ 

2 5x,'Cycle Time = ',f 12.3,' sec'/ 

3 5x,’%Onrime =',fl2.3,’%’/ 

4 5x,' Vented Mass = ',£12.3,' Ibm'/ 

5 5x, 'Specific Fuel Consumption = ',f 12.3,' Ibrn/hf) 
c 

c output listing 
c 

call prtout 
c 

c output plotting 
c 

if (opu .eq. 1.0) call crtplt(>l 1,12,2, 00, 0,0,subhd,xtitl, 

1 ypu,0.,0.,0.,0.,j,timep,pup,0,0) 

if (opl .eq. 1.0) call crtplt(12,12,2,00,0,0,subhd,xtitl, 

1 ypl,0.,0.,0.,0.j,timep,plp,0,0) 

if (otu .eq. 1.0) call crtplt(21, 00,2,00, 0,0,subhd,xtitl, 
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1 ytu,0.,0.,0.,0.j,timep,tup,0,0) 

if (otl .eq. 1.0) call crtplt(22, 00, 2,00,0,0, subhd.xtitl, 

1 ytl,0.,0.,0.,0.,j,timep,tlp,0,0) 

if (omu .eq. 1.0) call crtplt(-ll,22,2,00,0,0,subhd,xtitl, 

1 ymu,0.,0.,0.,0.j,timep,mup,0,0) 

if (oml .eq. 1.0) call crtplt( 12,00, 2,00,0,0,subhd,xtitl, 

1 yml,0.,0.,0.,0.j,timep,mlp,0,0) 

if (omwl .eq. 1.0) call crtplt(21,00,2,00,0,0,subhd,xtitl, 

1 ymwl,0.,0.,0.,0.j,timep,mwlp,0,0) 

if (omdsu .eq. 1.0) call crtplt(-ll,22,2,00,0,0,subhd,xtitl, 

1 ymdsu,0.,0.,0.,0.,j,timep,mdsup,0,0) 

if (omddu .eq. 1.0) call crtplt(12,00,2,00,0,0,subhd,xtitl, 

1 ymddu,0.,0.,0.,0.,j,timep,mddup,0,0) 

if (omdbw .eq. 1.0) call crtplt(21,00,2,00,0,0,subhd,xtitl, 

1 ymdbw,0.,0.,0.,0.j,timep,mdbwp,0,0) 

if (omdlu .eq. 1.0) call crtplt(22,00,2,10000,jsymb(l),0,subhd, 
1 xtitl,ymdif,0.,0.,0. ,0. j,timep,mdlup,0,0) 

if (omdul .eq. 1.0) call crtplt(00,00,0,10001 jsymb(2),0,0,0, 

1 0,0.,0.,0.,0.j,timep,mdulp,0,0) 

if (omdlu .eq. 1.0 .and. omdul .eq. 1.0) 

1 call crtkey(2jsymb,ketitl,- 1,-1) 

if (omds .eq. 1.0) call crtplt(-ll,12,2,00,0,0,subhd,xtitl, 

1 ymds,0.,0.,0.,0.,j,timep,mdsp,0,0) 

if (omdv .eq. 1.0) call crq>lt(12,12,2,00,0,0,subhd,xtitl, 

1 ymdv, 0 ., 0 ., 0 ., 0 .j,timep 4 ndvp, 0 , 0 ) 

if (onpump .eq. 1.0) call crtplt(-ll,12,2,00,0,0,subhd,xtiti, 

1 ynpump,0.,0.,0.,0.j,timep,npumpp,0,0) 

if (otw .eq. 1.0) call citplt(21,00,2,00,0,0,subhd,xtitl, 

1 ytw,0.,0.,0.,0.,j,timep,twp,0,0) 

if (odppmp .eq. 1.0) call crtplt(12,12iz,00,0,0,subhd,xtitl, 

1 ydppmp,0.,0.,0.,0.j,timcp,dppmp,0,0) 

stop 

end 

subroutine aiea(vl,y,awu,aul,awl,zdm) 
c 

c subroutine arca.f to calculate the heat transfer areas, 
c and liquid and gas heights of an elliptical bulkhead tank 
c 

reall 

common /tankdini/do,l,ho,t,cl 
common /tankvol/vl,v2,vbi,vbo,vci,vco 
common Aankarea/ab,ac, at 
data tol,nlim/.01,40/ 
r = do/2. - 1 
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h = ho - 1 
vb = vbi 
VC = vci 


c 

c liquid level is in upper bulkhead 
c 

if (vl .gt, v2) then 
V = vl - vb - VC 
do i = l,nlim 

fy = y*y*y - 3.*h*h*y + 3.*v*h*h/(3.14159*r*r) 
dfy = 3.*y*y - 3.*h*h 
dely = fy/dfy 
y = y - dely 

if (abs(dely) .le. tol) go to 5 
enddo 

5 z=y+h+l 
awl = ab + ac 

1 + 3.14159*(y*(cl*y*y + r*r)**.5 + r*r/cl** 5 

2 *(log(y*cl**.5 + (cl*y*y + r*r)**.5) - log(r))) 
ml = r*(l, - y/h) 

aul = 3.14159*ral*rul 
c 

c liquid level is in cylindrical segment 
c 

else 

if (vl .gt vl) then 
z = (vl - vb)/(3.14159*r*r) + h 
awl = ab + 2.*3.14159*r*(z - h) 
ml «r 

aul = 3.14159*rul*ml 
c 

c liquid level is in lower bulkhead 
c 

else 

if (vl .gt 0.) then 
V = vb - vl 
do i = l^Uim 

fy s y*y*y - 3.*h*h*y + 3.*v*h*h/(3.14159*r*r) 
dfy = 3.*y*y-3.*h*h 
dely = fy/dfy 
y = y - dely 

if (abs(dely) .le. tol) go to 10 
enddo 

10 z = h-y 
else 


m-73 



z = 0. 
endif 

awl = 3.14159*(h*(cl*h*h + r*r)**.5 

1 - (h - z)*(cl*(h - z)*(h - z) + r*r)**.5 

2 + r*r/cl**.5*(log(h*cl**.5 + (cl*h*h + r*r)**.5) 

3 - log((h - z)*cl**.5 + (cl*(h - z)*(h - z) + r*r)**.5))) 
ml = z*r/h 

aul = 3.14159*rul*ml 
endif 
endif 

awu = at - awl 
dul = 2.*nil 
hu = 2.*h + 1 - z 
return 
end 

subroutine frict(d,ed,mdot,visc,f) 
c 

c subroutine Mctf to calculate the friction coefficient 
c for flow in a pipe 
c 

real mdot 

re = 4,*mdot/(3.14159*visc*d) 
if (ic .It 2300.) f = 64./(re + 1.) 
if (re .gt 2300.) 

1 f « 0.25/Gogl0(ed/3.7 + 2.51/(rc*sqrt(.0056 + 

2 .5/(re**.32)))))**2 
return 

end 

subroutine htc(tl,t2,d,cp^o,beta4t,mu,a,h) 
c 

c This subroutine computes the free convection heat-transfer 
c coefficient for horizontal and vertical surfaces 
c 

real ltmu,nu 
if (d .eq. 0.) h = 0. 
if (d .eq. 0.) return 

ra = 3600*32.2*a*beta*abs(tl - t2)*d**3*rho*rho*cp 
1 /(mu*k) 

nu = 0.555*ra**0.25 + 0.447 

h = nu*k/d 

return 

end 
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subroutine pdrop(mdot,rho,mu,U 41,d2,dptot) 
c 

c subroutine pdrop.f to calculate the pressure drops 
c between the pump outlet and spray manifold inlet 
c 

real reno(9)ddossl(9)4doss2(9) 

real mdot,mudldcbldcb2^b34cb4,kcJcflm 

data reno/1 .5e+5,2.0e+5,3*0e+5,4*0e+5,6‘0e+5,8.0e+5, l.Oe+6, 

1 2.0e+6,3*0e+6/ 

data klossi/0.32,0.26,0.21,0.19, 0.173, 0.168,0.163,0.160, 0.158/ 
data kloss2/0.20,0. 15,0. 14,0. 128,0. 12,0. 1 1 8,0. 1 17,0. 1 15,0. 1 14/ 
ed = l.Oe-6 
gc = 32.2 

al = 3.14159*dl*dl/4.0 
a2 = 3.14159*d2*d2/4.0 
c 

c 90-degree bend at pump outlet 
c 

rel = 4.0*mdot/(3.14159*mu*dl) 
call value(9,reno4clossl,re 1 Jcb 1) 
c 

c straight section downstream of 90-degree bend 
c 

call ffict(dl,ed,mdot,mu,fl) 
kl=fl*ll/dl 
c 

c reducer 
c 

kc = 0.5*(1.0-(d2/dl)**2) 
c 

c 132.5-degree bend 
c 

re2 = 4.0*mdot/(3. 14159*mu*d2) 
call value(9,renoJdossl,re2Jcb2) 
kb2=1.22*kb2 
c 

c flowmeter 
c 

kflm= 1.308 
c 

c 95.5-degree bend downstream of flowmeter 
c 

call value(9,reno4doss24e24cb3) 
kb3 = 1.03*kb3 
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c 48-degree bend 
c 

call value(94«no,kloss24:e24cb4) 
kb4 = 0.66*kb4 
c 

c pressure drops 
c 

const = mdot*mdot/(2.0*rho*gc* 144.0) 

dpbl = kbl*const/al**2 

dpi =kl*consl/al**2 

dpc = kc*const/a2**2 

dpb2 = kb2*const/a2**2 

dpflm = kflm*const/a2**2 

dpb3 = kb3*const/a2**2 

dpb4 = kb4*const/a2**2 

dptot = dpbl + dpi + dpc + dpb2 + dpflm + dpb3 + dpb4 
c write (6, 12) dpb 1 ,dp 1 ,dpc,dpb2,dpflm,dpb3,dpb4,dptot 
cl2 format (3x, 'dpbl =’48.5,3x,'dpl =',f8.5,3x,’dpc =‘^.5/ 
c 1 3x/dpb2 = ’^.53x,’dpflm = ’^.5,3x,'dpb3 = *^.5/ 

c 2 3x,'dpb4 = ’48.5^x, ’dptot = ',f8.5) 

return 
end 

subroutine prtout 
c 

c subroutine prtoutf to print the complete output of 
c program tvs.f in a file named outdat 
c 

real mu^nl^nwl^nv 

real mdlu,mds 4 ndsl 4 ndsu 4 nddu 4 ndbw,mdul,mdcond 4 ipump 
conunon /tankout/nstep4iline,iprint,time(3000), 

1 pu(3000),tu(3000),vu(3000),mu(3000),tw(3000), 

2 pl(3000),tl(3000).vl(3000),ml(3000), 

3 pwl(3000),twl(3000),vwl(3000),mwl(3000), 

4 mv(3000),ts(3000)4ndlu(3000),mds(3000),mdsl(3000), 

5 mdsu(3000),mddu(3000),mdbw(3000)4ndul(3000),mdcond(3000), 

6 qwu(3000),quwl(3000),qwl(3000),qul(3000),qud(3000), 

7 qus(3000),qls(3000),npump(3000),dppump(3000) 
open (unit=15,file=5'outdat') 

icount - nstep/(nline*iprint) + 1 
do j = Idcount 
imin = imax*f 1 
imax = j*nline*iprint 
if (imax .gt nstep) imax = nstep 
write (15,1) 


Iff 

E 

I 


i 
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1 format (l\t6;Time\tl8,'pU',t28;TU*,t38,W,t48;MU’, 

1 t58;TW’,t68;pL’,t78;TL',t88;VL',t98/ML’, 

2 tl07,'pWL',tl 17,'TWL\tl27/VWL’,tl37;MWL’, 

3 tl48,’MV',tl58,TS’,tl66;mdLU’) 
write (15,2) 

2 format (t7,'sec',tl6,'psia',t29,'R’,t37,'ft3',t47,'lbm’, 

1 t59,‘R\t66,‘psia',t79,'R',t87,'ft3’,t97,'lbm’, 

2 tl06,'psia',tl 19,‘R',tl27,’ft3',tl37,‘lbm’, 

3 tl47,'lbm',tl59;R',tl63,lbm/sec'/) 
do 20 i = imin4max 

if (mod(i-l,iprint) .ne. 0.0) go to 20 
write (15,3) time(i),pu(i),tu(i),vu(i),mu(i), 

1 tw(i),pl(i),tl(i),vl(i),ml(i), 

2 pwl(i),twl(i),vwl(i),mwl(i), 

3 mv(i),ts(i),mdlu(i) 

3 format (f9.1,15(lx,f9.3),lx,c9.3) 
c3 . format (fl0.1,15(lx,el0.4)) 

20 continue 

write (15,4) 

4 format Cl’,t6,Time',tl7,’mdS’,t26;mdSL’,t36,’mdSU’, 

1 t46,’mdDir,t56;mdBW,t66,’mdUL‘,t74,’mdCOND‘, 

2 t87;qW,t96,’qUWL’,tl07,'qWL‘,tl 17,’qUL’, 

3 tl27,’qUD',tl37,’qUS’,tl47,’qLS', 

4 tl55,'Npump',tl64,'dppump') 
write (15,5) 

5 format (t7,'sec',tl3,’lbm/sec',t23,’lbm/sec',t33,lbm/sec', 

1 t43,'lbm/sec',t53,lbm/sec',t63,'lbm/sec',t73,lbm/sec', 

2 t84,'Btu/hri,t94;Btu/hr‘,tl04,'Btu/h^ 14,’Btu/hr', 

3 tl24;Btu/hr',tl34,*Btu/hr’,tl44,'BtiVhr\ 

4 tl57,'rpm',tl67,’psi’/) 
do 30 i = imin4max 

if (mod(i-l,iprint) .ne. 0.0) go to 30 
write (15,6) time(i),mds(i),mdsl(i),mdsu(i), 

1 mddu(i) 4 ndbw(i),mdul(i),mdcond(i), 

2 qwu(i),quwl(i),qwl(i),qul(i), 

3 qud(i),qus(i),qls(i), 

4 npump(i),dppump(i) 

6 format (f9.1,7(lx,c9.3),9(lx,f9.3)) 

30 continue 

enddo 

return 

end 

subroutine spray(pman,mdoti,pull,zu,zliq,ztank,dchar,rhol,acc,ed, 
1 time,prtsp,mdot,mdsu,mdsl,veld,ndrop,norfu) 
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c 

c subroutine spray .f to model flow in the spray injection tube 
c 

real rod(12),lod(12),roverd,loverd 

real mdin(l 00),mdout( 100),mds( 1 00),as( 1 00),asp( 1 00) 

real vd(100),veld(100),ndrop(100) 

real pin( 1 00),pout( 1 00),pnode( 1 00),ptank( 1 00),x(l 00),delpt(200) 

realmdotimdoti,mdotsi,mdoto*mdots,mdsu,mdsl 

real kfsnvkbsm^csnU^bsiJ^ 

common /sprayin/dsm,zsm,dsi,zsi,n,nbar, 

1 ks,cds,as,vd/overd,dmdt,tol,nlim, 

2 rho,visc 

data nb/12/rod/l., 1.5, 2.,3.,4.,6.,8.,10.,12.,14.,16.,20./ 

datalod/20.,14.,12.,12.,14.,17.,24.,30„34.,38.,42.,50./ 

psmi = 144.*pman 

pu = 144.*pull 

dmdot = dmdt 

dzsi = zsi/n 

asm = 3.14159*dsm*dsm/4. 
asi = 3. 14159*dsi*dsi/4. 
kcsm = .5*(1. - (dsi/dsm)**2) 
call value(nb,rod,lod,ioverd4overd) 
gc = 32.2 
mdot = mdoti 
do 10 j = l,nlim 

qsm = (mdot/asm)**2/(2.*rho*gc) 
call &ict(dsm,ed,mdot,visc,fsm) 
kfsm s fsm*zsm/dsm 
kbsm = fsm*loverd 
dpfsm = qsm*(kfsm + kbsm + kcsm) 
phsm ss rho*acc*zsm 
dpsm = dpfsm + phsm 
psmo = psmi • dpsm 
mdotsi = mdot/nbar 
qi = (mdotsi/asi)**2/(2.*rho*gc) 
call frict(dsi,ecMndotsi,visc/si) 
kbsi = fsi*loverd 
pi = psmo-qi*kbsi 
mdsu - 0. 
mdsl = 0. 
norfu = 0 
do i s l,n 

x(i) = i*dzsi - dzsi/2. 
aorf = as(i) 
z = ztank - zsi + x(i) 
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if (z .It zu) pt = pu 

if (z .ge. zu) pt = pu + rhol*acc*(z - zu) 
callpres(i,n,pi,po,pn,pt,dpf»ph,nidotsi,mdoto,mdots, 

1 dzsi,dsi,ed,asi,aorf,ks,rho,visc,acc) 

pin(i)=pi 
pout(i) = po 
pnode(i) = pn 
ptank(i) = pt 
mdin(i) = mdotsi 
mdout(i) = mdoto 
pi = pout(i) 
mdotsi = mdout(i) 
mds(i) = mdots 

if (aoif .gt 0.0) veld(i) = mdots/(rho’*‘cds*aorf) 
if (aoif .le. 0.0) veld(i) = 0.0 
if (veld(i) .gt 0.0) 

1 ndrop(i) = nbar*mdots*dchar/(rho*vd(i)*veld(i)) 
if (veld(i) .le. 0.0) ndrop(i) = 0.0 
if (z .It zu) noifu = norfu + 1 
if (z .It zu) mdsu = mdsu + mdots 
if (z .ge. zu) mdsl = mdsl + mdots 
enddo 

mdsu = nbar*mdsu 
mdsl = nbar*mdsl 

ptcal = pn - ks/(2.*rho*gc)*(mdots/aorf)**2 

delpt(j) = ptcal - pt 

dpt =s delpt(j)/144.0 

if (abs(dpt) .It tol) go to 15 

if (j .eq. 1) then 

if (delpt(j) .It 0.) mdot - mdot - dmdot 
if (delpt(j) .gt 0.) mdot = mdot + dmdot 
else 

prod s= delpt(i)*delpt(j- 1) 

if (delpt(j) .It 0. .and. prod .gt 0.) mdot « mdot - dmdot 
if (delpt(j) .gt 0. .and. prod .gt 0.) mdot = mdot + dmdot 
if (delpt(j) .It 0. .and. prod .It 0.) then 
dmdot = dmdot/2. 
mdot = mdot - dmdot 
endif 

if (delpt(j) .gt 0. .and. prod .It 0.) then 
dmdot = dmdot/2. 
mdot = mdot + dmdot 
endif 
endif 

10 continue 



if (abs(dpt) .gt tol) write (6,1002) time,dpt 
1002 format C*** tank pressure does not converge at time * ',fl0.2, 

1 ' sec,’,' delpt = ',e9.3,' psi ***’) 

15 if (time .ge. prtsp .and. time .le. (prtsp + 0. 1)) then 
. dppump = (psmi - pt)/144. 
dpsm = dpsm/144. 
dpsi = (pi - pt)/144. 
hliq =s 12.*zliq 
write (6,1) 

1 format (/5x,’SPRAY MANIFOLD/INJECnON TUBE FLOW MODEL'/) 
write (6,2) acc,puU,pman4iliq,rho,visc 

2 format (5x,‘ Acceleration Level ',f6.1,’g'/ 

1 5x,'UUage Pressure *,f6.3,' psia’/ 

2 5x,'Spray Manifold Inlet Pressure ',f6.3,' psia'/ 

3 5x,'Liquid Level ',f6.1,’in'/ 

4 5x,'Liquid Density ',f6.3,' lbm/ft3’/ 

5 5x,‘Liquid Viscosity ',e9.3,' Ibm/ft-sec*/) 

write (6,3) mdot,dppump,dpsm,dpsi 

3 format (5x,’Pump Flow Rate ',f6.3,' Ibm/sec'/ 

1 5x,'Pump Pressure Rise ',f6.3,' psi*/ 

2 5x,’Spray Manifold Tube Delta p '4^6.3,' psi'/ 

3 5x,'Spray Injection Tube Delta p '46.3,' psi'/) 
write (6,4) 

4 format (5x,'Node',2x,*Distance',3x,Tnlef ,2x,'Outlet*, 

1 3x,'Nodal',3x,Tank',5x,'Inlet',4x,'Outlet', 

2 3x,'Injection',3x,'Injection*,3x,'Orifice’) 
write (6,5) 

5 format (24x,'p',7x,'p’,7x,'p',7x,'p’, 

1 7x,'mdot'4x,*mdot',7x,'mdof,6x,*Velocity', 

2 5x,'CdA') 
write (6,6) 

6 format (13x,*(in)',4x,’(psia)',2x,’(psia)',2x,'(psia)’,2x,'(psia)’, 

1 2x,'(lbm/scc)',lx,'(lbin/sec)',lx,'Gbm/sec)'4x,'(^s)', 

2 6x,'(in2)'/) 
do i s l,n 

pin(i) = pin(i)/144. 
pout(i) = pout(i)/144. 
pnode(i) = pnode(i)/144. 
ptank(i) = ptank(i)/144. 
asp(i) « 144.*as(i) 
x(i) = 12.*x(i) 

write (6,7) i,x(i),pin(i),pout(i),pnode(i),ptank(i), 

1 mdin(i) 4 ndout® 4 nds(i),veld(i),asp(i) 

7 format (5x,i3,4x46.2,3x,4(f6.3,2x), 

1 3(e9.3,lx),4x,f5.2,4x48.6) 



enddo 

endif 

return 

end 

subroutine value(np,x,y,xin,yout) 
c 

c This subroutine performs lagrangian interpolation within a set 
c of (x,y) pairs to give yout corresponding to xin 
c 

dimension x(np),y(np) 
if (xin .le. x(l)) yout = y(l) 
if (xin .le. x(l)) return 
if (xin .ge. x(np)) yout = y(np) 
if (xin .ge. x(np)) return 
do 10 i = l,np 
k = i 

if (xin .It x(i)) go to 30 
10 continue 

30 ffr = (xin-x(k-l))/(x(k)-x(k-l)) 
yout = y(k - 1) + ffr*(y(k) - y(k -1)) 
return 
end 

SUBROUTINE VALUE2(NPX,NPY,X,Y^,XIN,YIN^OUT) 
DIMENSION X(NPX),Y(NPY),Z(NPX,NPY) 

C 

IF(XIN .LE. X(l) .AND. YIN .LE. Y(l))ZOUT=Z(l,l) 
worn .LE. X(l) .AND. YIN .LE. Y(1))RETURN 
C 

IF(XIN .GE. X(NPX) .AND. YIN .LE. Y(l))ZOUT=Z(NPX,l) 

IF(XIN .GE. X(NPX) .AND. YIN .LE. Y(1))RETURN 
C 

IF(XIN .LE. X(l) .AND. YIN .GE. Y(NPY))ZOUT=Z(l,NPY) 

IF(XIN .LE. X(l) .AND. YIN .GE. Y(NPY))RETURN 
C 

IF(XIN .GE. X(NPX) .AND. YIN .GE. Y(NPY))ZOUT=Z(NPX,NPY) 
IF(XIN .GE. X(NPX) .AND. YIN .GE. Y(NPY))RETURN 
C 

IF(XIN .GT. X(l))GO TO 30 
DO20I=lJ^PY 
M=I 

IF(YIN .LT. Y(I))GO TO 25 
20 CONTINUE 

25 FFRY=(YIN-Y(M-1))/(Y(M)-Y(M-1)) 
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r 

I 


Z0UT=Z(1,M-1)+FFRY*(Z(1,M)-Z(1,M-1)) 

RETURN 

C 

30 IF(XIN .LT. X(NPX))GO TO 60 
DO 50 1=1, NPY 
M=I 

IF(YIN .LT. Y(I))GO TO 55 
50 CONTINUE 

55 FFRY=(YIN-Y(M-1))/(Y(M)-Y(M-1)) 
ZOUT=Z(NPX,M-l)+FFRY*(Z(NPX,M)-Z(NPX>I-l)) 
RETURN 
C 

60 IF(YIN .GT. Y(l))GO TO 90 
DO 80 I=1,NPX 
L=I 

IF(XIN .LT. X(I))GO TO 85 
80 CONTINUE 

85 FFRX=(XIN-X(L-1))/(X(L)-X(L-1)) 
ZOUT=Z(L-l,l)+FFRX*(Z(L,l)-Z(L-U)) 

RETURN 

C 

90 IF(YIN .LT. Y(NPY))GO TO 120 
DO110I=U^PX 
L=I 

W(xm .LT. X(I))GO TO 115 
110 CONTINUE 

115 FFRX=(XIN-X(L-1))/(X(L)-X(L-1)) 
ZOUT=Z(L-l,NPY)+FFRX*(Z(L,NPY)-Z(L-l,NPY)) 
RETURN 
C 

120D0130I=WX 

L=I 

IF(XIN .LT. X(I))GO TO 135 
130 CONTINUE 

135 FXR=(XIN-X(L-1))/(X(L)-X(L-1)) 

DO 140I*1,NPY 
M=I 

IF(YIN .LT. Y(D)GO TO 145 
140 CONTINUE 

145 FYR=(YIN-Y(M-l))/Cy(M)-Y(M-l)) 

C 

ZXLO=Z(L-l^-l)+FYR*(Z(L-l>I)-Z(L-lM-l)) 

ZXHI=Z(L>I-1)+FYR*(Z(L,M)-Z(L>I-1)) 

C 

ZOUT=ZXLO+FXR*(ZXHI-ZXLO) 


i 
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RETURN 

END 


subroutine volarea(vt) 

subroutine volarea.f to calculate the volumes and areas of 
an elliptical bulkhead tank 

reall 

common /tankdim/do,l,ho,t,cl 

common /tankvol/vl,v2,vbi,vbo,vci,vco 

common /tankaiea/abi,aci,ati 

ro = do/2. 

ri = ro - 1 

hi = ho - 1 

xkl =ri/hi 

cl = xkl**4 - xkl*xkl 

internal bulkhead 

vbi = 2./3.*3.14159*ri*ri*hi 

abi = 3.14159*(hi*(cl*hi*hi + ri*ri)** 5 

1 + ri*ri/cl**.5*(log(W*cl**.5 

2 + (cl*hi*hi + ri*ri)**.5) - log(ri))) 

external bulkhead 

vbo = 2./3.*3.14159*ro*ro*ho 

abo = 3.14159*(ho*(cl*ho*ho + ro*ro)**.5 

1 + ro*ro/cl**.5*(log(ho*cl**.5 

2 + (cl *ho*ho + io*ro)**.5) - log(ro))) 

internal cylinder 

vci = 3. 14159*ii*ii*l 
aci = 2.*3.14159*ri*l 

external cylinder 

vco = 3.14159*ro*ro*l 
aco = 2.*3.14159*ro*l 

tank 

vl = vbi 
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v2 = vbi + vci 
vt = 2.*vbi + vci 
ati = 2.*abi + aci 
return 
end 

subroutine vwall(hliq,vtw) 
c 

c subroutine vwall.f to calculate the wall volume exposed to 
c the ullage gas 
c 

reall 

common /tankdim/do,l,ho,t,cl 
common /tankvolA^l,v2,vbi,vbo,vci,vco 
ro = do/2, 
ri = ro - 1 
hi = ho - 1 
vbw = vbo - vbi 
vcw = vco - vci 
if (hliq .gt (hi + 1)) then 
h = hliq-(hi + l) 

vi = 3.14159*h*(ri/hi)**2*(hi*hi - h*h/3.) 

VO = 3.14159*h*(ro/ho)**2*(ho*ho - h*h/3.) 

V = VO - vi 
vtw = vbw - v 
else 

if (hliq .gt hi) then 
h - hi + 1 - hliq 
vi = 3.14159*ri*ii*h 
vo = 3.14159*ro*ro*h 

V = VO - vi 
vtw = vbw + v 

else 

h = hi ’ hliq 

vi 3.14159*h*(ri/hi)**2*(hi*hi - h*h/3.) 

VO = 3.14159*h*(ro/ho)**2*(ho*ho - h*h/3.) 

V = vo - vi 

vtw * vbw + vcw + V 
endif 
endif 
return 
end 




m 
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3.5 Input Data 

3.5.1 Heat Exchanger Model 


MASSICI PIC I PSTI IMDOT(VNT)l PBACK I Program is vent2 

* % I ♦ ♦ ♦ 9k 4c ♦ 3(( % I * >i( >l« 9k ♦ 3k 9k 1 4( 9|( 9k 4c)k 3i( 1 3|« )k * * ♦ 3k 4* * I 

270.0 20.590 20.645 0.00475 19.5 

PROP I PTP I AVLIQ I PAMB I G I 

1.0 1.021 3500.0 0.0 0.0 

M I VT I AX I VA I Twall I 
25 68.9342 0.04500 0 32.00 

MDOT(sp)IDI(sp)ITHKNESSIVFLWDII E/D I 

4iak9k9k9kik3k9k9k|3k9k9k9k9k3k9k3k9k|3k3kik9k9k9k3k9k9k|9k9k9k9k9k9k3k3k9k|3k9k9k9k3k9k9k9k3k| 

,0.32 1.18 0.035 0.25 l.Oe-6 

9k 9k 9k 9k 9k 9k )k 3k 9k 1 9k 9k 9k 9k 9k 9k 9k 9k 9k 1 9k 9k 9k 3k 9k 9k 9k 9k 3k 1 3k 3k 9k 9k 9k 9k 3k ik 9k I )k 9k 9k 9k 9k 9k 9k 3k 9k I 

DPINC I ERRMX I PERRMXI DEBUG lEQDIAMI 

9k 9k 9k 9k 9k 9k 9k 9k 9k 1 9k 9k 9k 9k 3k 9k 9k 9k 9k 1 9k 9k 9k 9k 9k 9k 3k 9k 9k 1 9k 9k 9k 9k 3k 9k 9k 9k 9k 1 9k 9k 9k 9k 9k 9k 3k 9k 9k I 

0.005 0.005 0.0010 0.0 0.134 

9k9k9k9k9k9k9k9k9k|9k9k9k9k9k9k9k9k9k|9k9k9k3k9k>k9k9k9k|9k9k9k9k9k9k9k9k9k|9k9k9k>k9k9k9k9k9k| 

FINTIM I PRDEL I QDTERR I DELT I DELPTP I 

3k9k9k9k9k9k9k9k9k|9k9k9k9k9k9k9k9k9k|9k9k9k9k9k>k9kik9k}9k9k)k9k9k9k9k9k9k]9k9k9k9k9k9k9k9k9k| 

2.2 0.02 0.03 0.01 0.08 

9k 9k 9k 9k 9k 9k 9k 9k 9k 1 9k 3k 9k 9k 9k 9k 9k 9k 9k 1 9k 9k 9k 9k 9k 9k 9k 9k 9k 1 9k 9k 9k 9k 9k 3k 3k 9k 9k 1 9k 9k 3k 9k 9k 9k 9k 9k 9k I 

OT I OM I OP IDELQdotl I 

9k9k9k9k9k9k3k9k9k|3k9k9k9k9k9k9k3k9k|9k9k9k9k9k3k9k9k9k|9k9k3k9k9k9k3k9k9k]9k9k9k3k9k9k9k9k9k| 

0 0 0 0.0010 

*9k9k9k9k9k9k9k9kj9k9k9k9k9k3k9k9k9k|3k9k9k9k9k9k9k9k3k[9k9k9k9k9k9k9k9k9kj 

A(l-M) 1 DH(l-M) IQDOT(l-M)l LENGTH I 

9k3k9k9k9k3k9k9k9k]9k9k9k9k9k9k9k9k9kj3k9k9k9k9k9k9k9k9k{3k3k9k9k9k3k9k9k9kj 


0.04500 

0.0 

1.00 

0.0 

0.27721 

0.0 

1.00 

6.0 

0.27721 

0.0 

0.975 

6.0 

0.27721 

0.0 

0.95 

6.0 

0.27721 

0.0 

0.925 

6.0 

0.27721 

0.0 

0.90 

6.0 

0.27721 

0.0 

0.875 

6.0 

0.27721 

0.0 

0.85 

6.0 

0.27721 

0.0 

0.825 

6.0 

0.27721 

0.0 

0.80 

6.0 


0.27721 

0.0 

0.75 

6.0 

0.27721 

0.0 

0.70 

6.0 

0.27721 

0.0 

0.65 

6.0 

0.27721 

0.0 

0.60 

6.0 

0.27721 

0.0 

0.500 

6.0 

0.27721 

0.0 

0.450 

4.0 

0.27721 

0.0 

0.400 

4.0 

0.27721 

0.0 

0.350 

4.0 

0.27721 

0.0 

0.300 

4.0 

0.27721 

0.0 

0.250 

4.0 

0.27721 

0.0 

0.200 

4.0 

0.27721 

0.0 

0.150 

4.0 

0.27721 

0.0 

0.100 

4.0 

0.27721 

0.0 

0.050 

4.0 

0.00372 

0.0 

0.000 

0.0 


Simulation of LH2 Vent Thru Zero g Vent System Heat Exchanger (4/6/93) 
Baseline Vent Area of 0.00372 in2 (Mdot & Qdot Trade) delta q = 0.0010 

H2 MISC DATA 
LMXPAR 
18 


-2.0 

2.107 

-1.699 

1.835 

-1.3979 

1.585 

-1.1549 

1.3874 

-1.0 

1.2672 

-0.6990 

1.0492 

-0.3979 

0.8482 

-0.1549 

0.7024 

0.0 

0.6232 

0.3010 

0.4914 

0.6020 

0.3766 

0.8451 

0.2923 

1.0 

0.2430 

1.301 

0.1703 

1.602 

0.1106 

1.8451 

0.0682 

2.0 

0.0453 

10.0 

0.04 

LVISP 


31 


0.0 

1.985e-5 

1.021 

1.751e-5 
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2.553 

1.4e-5 

4.17 

1.24e-5 

6.446 

l.lle-5 

9.527 

1.006e-5 

13.561 

0.917e-5 

15.984 

0.877e-5 

18.694 

0.84e-5 

21.723 

0.805e-5 

25.089 

0.773e-5 

28.813 

0.742C-5 

32.915 

0.713e-5 

37.415 

0.685e-5 

42.334 

0.659C-5 

47.693 

0.634C-5 

53.514 

0.609e-5 

59.817 

0.586C-5 

66.625 

0.563C-5 

73.957 

0.54C-5 

81.838 

0.519C-5 

90.287 

0.497e-5 

99.329 

0.476C-5 

108.987 

0.454C-5 

119.297 

0.433C-5 

130.299 

0.410e-5 

142.027 

0.387C-5 

154.522 

0.362C-5 

167.848 

0.333e-5 

182.136 

0.292e-5 

187.510 

0.24e-5 

GVISP 


24 


0.0 

0.045e-5 

1.021 

0.05e-5 

2.553 

0.057C-5 

9.527 

0.07C-5 

18.694 

0.079C-5 

32.915 

0.089C-5 

37.415 

0.092C-5 

42.334 

0.095C-5 

47.693 

0.097C-5 

53.514 

O.le-5 

59.817 

0.103e-5 

66.625 

0.106C-5 

73.957 

0.109e-5 



81.838 

0.112e-5 


90.287 

0.116C-5 

f" 

99.329 

0.12e-5 

1- 

(■v- 

108.987 

0.124C-5 


119.297 

0.129e-5 

f 

130.299 

0.135C-5 

1 

142.027 

0.143e-5 


154.522 

0.153e-5 

p 

167.848 

0.168e-5 

I 

182.136 

0.196e-5 


187.51 

0.24C-5 

&■ 

f-''' 


H2 RHO DATA 

fee 

LRHOP 


i 

39 



0.0325 

5.385 


1.01 

5.385 


1.021 

4.80827 

1: 

1.073 

4.80377 


1.462 

4.77434 


1.950 

4.74423 


2.553 

4.71359 


3.288 

4.68219 


4.170 

4.65003 

k. 

5.217 

4.61705 


6.446 

4.58321 

L\ 

7.877 

4.54844 

tiv. 

9.527 

4.51266 


11.416 

4.47582 

5 :' 

13.561 

4.43782 


15.984 

4.39858 


18.694 

4.35801 


21.723 

4.31600 


25.089 

4.27243 

, . 

28.813 

4.22718 

k 

32.915 

4.18010 


37.415 

4.13102 


42.334 

4.07975 


47.693 

4.02608 


53.514 

3.96975 


59.817 

3.9105 


66.625 

3.8479 


73.957 

3.7815 

1: 

81.838 

3,7108 


90.287 

3.635 
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99.33 

3.5534 

108.987 

3.4646 

119.297 

3.3668 

130.299 

3.2572 

142.027 

3.1314 

154.522 

2.9808 

167.848 

2.7851 

182.136 

2.4567 

187.51 

1.9620 

GRHOP 


40 


0.0325 

0.000345 

0.1137 

0.00108 

0.3163 

0.002742 

1.021 

0.00784 

1.073 

0.00819 

1.462 

0.01077 

1.950 

0.01391 

2.553 

0.01765 

3.288 

0.02207 

4.170 

0.02724 

5.217 

0.03322 

6.446 

0.04008 

7.877 

0.04790 

9.527 

0.05675 

11.416 

0.06671 

13.561 

0.07787 

15.984 

0.09011 

18.694 

0.10395 

21.723 

0.11930 

25.089 

0.13629 

28.813 

0.15504 

32.915 

0.17568 

37.415 

0.19839 

42.334 

0.22335 

47.693 

0.25078 

53.514 

0.28094 

59.817 

0.3141 

66.625 

0.35076 

73.957 

0.39126 

81.838 

0.43621 

90.287 

0.48635 

99.329 

0.54266 

108.987 

0.60647 
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119.297 

0.67967 

130.299 

0.76511 

142.027 

0.86751 

154.522 

0.99566 

167.848 

1.17023 

182.136 

1.4779 

187.510 

1.96202 

RHOVP 

20 

0.000345 

0.0325 

0.00108 

0.1137 

0.002742 

0.3163 

0.00784 

1.021 

0.00819 

1.073 

0.01077 

1.462 

0.01391 

1.950 

0,01765 

2.553 

0.02207 

3.288 

0.02724 

4.170 

0.03322 

5.217 

0.04008 

6.446 

0.0479 

7.877 

0.05675 

9.527 

0.06671 

11.416 

0.07787 

13.561 

0.09011 

15.984 

0.10395 

18.694 

0.11930 

21.723 

0.13629 

25.089 

SCLENT 
10 7 

10.0 3 

34.0 

-115.878 

34.263 

-115.318 

36.483 

-110.59 

14.7 4 

34.0 

-115.740 

36.0 

-111.351 

36.483 

-110.240 

36.608 

-109.95 

15.0 4 

34.0 

-115.731 

36.0 

-111.343 
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36.608 

- 109.942 

38.444 

- 105.71 

20.0 

5 

34.0 

- 115.583 

36.0 

- 111.201 

38.0 

- 106.489 

38.444 

- 105.395 

39.975 

- 101.62 

25.0 

5 

34.0 

- 115.435 

36.0 

- 111.059 

38.0 

- 106.356 

39.975 

- 101.358 

41.299 

- 98.01 

30.0 

6 

34.0 

- 115.287 

36.0 

- 110.917 

38.0 

- 106.223 

40.0 

- 101.170 

41.299 

- 97.667 

42.475 

- 94.5 

35.0 

7 

34.0 

- 115.139 

36.0 

- 110.775 

38.0 

- 106.090 

40.0 

- 101.052 

42.0 

- 95.584 

42.475 

- 94.227 

43.536 

- 91.2 

40.0 

7 

34.0 

- 114.991 

36.0 

- 110.633 

38.0 

- 105.957 

40.0 

- 100.929 

42.0 

- 95.477 

43.536 

- 90.977 

45.406 

- 85.50 

45.0 

7 

34.0 

- 114.842 

36.0 

- 110.491 

38.0 

- 105.823 

40.0 

- 100.806 

42.0 

- 95.370 

44.0 

- 89.454 

45.406 

- 85.08 
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50.0 

7 

34.0 

-114.694 

36.0 

-110.348 

38.0 

-105.689 

40.0 

-100.682 

42.0 

-95.262 

44.0 

-89.366 

45.406 

-84.890 


LH2 SATURATED PROPERTIES 


* Pressure 
16 

* Cp 

* Thermal Cond* Surface T * Sat Temp 

1.021 

1.557 

0,04199 

1.7076e-5 

24.845 

1.462 

1.619 

0.04551 

1.6469e-5 

26.0 

2.553 

1.723 

0.05004 

1.5419e-5 

28.0 

4.17 

1.849 

0.05297 

1.4374e-5 

30.0 

6.446 

1.985 

0.05489 

1.3333e-5 

32.0 

9.527 

2.125 

0.05601 

1.2298e-5 

34.0 

13.561 

2.27 

0.05690 

1.1267e-5 

36.0 

18.694 

2.443 

0.05793 

1.0243e-5 

38.0 

25.089 

2.637 

0.05843 

0.9225e-5 

40.0 

28.813 

2.743 

0.05851 

0.8718e-5 

41.0 

32.915 

2.848 

0.05848 

0.8213e-5 

42.0 

42.334 

3.097 

0.05811 

0.7209e-5 

44,0 

53.514 

3.393 

0.05734 

0.6214e-5 

46.0 

66.625 

3.772 

0.05616 

0.5228e-5 

48.0 

81.838 

4.307 

0.05459 

0.4253e-5 

50.0 

99.329 

5.074 

0.05258 

0.3292C-5 

52.0 


GH2 PROPERTIES 

* Cp ♦ Thermal Cond* Viscosity * 
6 14 50.0 

1.0 


2.511 

7.21e-3 

5.0C-7 

2.506 

7.35C-3 

5.2C-7 

2.499 

7.63e-3 

5.6C-7 

2.493 

7.92e-3 

6.1e-7 

2.489 

8.24C-3 

6.5e-7 

2.486 

8.56C-3 

6.9C-7 

2.483 

8.96C-3 

7.3e-7 

2.481 

9.47e-3 

7.7e-7 

2.479 

9.98e-3 

8.1e-7 

2.478 

1.049e-2 

8.5e-7 

2.477 

l.lOe-2 

8.9e-7 

2.475 

1.151C-2 

9.3e-7 

2.474 

1.202C-2 

9.7e-7 
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1.252e-2 


l.Oe-6 


2.474 

5.0 
2.652 
2.636 
2.610 
2.589 
2.574 
2.561 
2.551 
2.543 
2.536 
2.524 
2.524 
2.52 
2.516 
2.512 

10.0 
2.790 
2.768 
2.726 
2.696 
2.671 
2.649 
2.633 
2.619 
2.607 
2.596 
2.587 
2.578 
2.571 
2.564 

15.0 

2.912 

2.886 

2.837 

2.800 

2.766 

2.740 

2.717 

2.698 

2.681 

2.666 

2.653 

2.641 

2.631 


6.33e>3 

7.815e-3 

8.707e-3 

8.879e-3 

9.174e-3 

9.550e-3 

9.933e-3 

1.032e-2 

1.070e-2 

1.109C-2 

1.148C-2 

1.186C-2 

1.224C-2 

1.262C-2 

7.04C-3 

7.892e-3 

9.425C-3 

9.714C-3 

l.OOle-2 

1.030C-2 

1.061e-2 

1.092C-2 

1.122C-2 

1.153C-2 

U83e-2 

1.214C-2 

1.245C-2 

1.276C-2 

9.84C-3 

9.95C-3 

1.017e-2 

1.040C-2 

1.064C-2 

1.088e-2 

1.113C-2 

1.138C-2 

1.163C-2 

1.188C-2 

1.213e-2 

1.239C-2 

1.264C-2 


6.3e-7 

6.43e-7 

6.72C-7 

7.03e-7 

7.35e-7 

7.73e-7 

8.05e-7 

8.32e-7 

8.56e-7 

8.87C-7 

9.18e-7 

9.49C-7 

9.79C-7 

l.Ole-6 

7.1e-7 

7.21C-7 

7.43e-7 

7.68C-7 

7.93C-7 

8.19C-7 

8.38e-7 

8.61C-7 

8.93e-7 

9.19e-7 

9.43C-7 

9.62C-7 

9.85C-7 

l.Ole-6 

7.6C-7 

7.72C-7 

7.94C-7 

8.16e-7 

8.37C-7 

8.59C-7 

8.78C-7 

8.94C-7 

9.13C-7 

9.34C-7 

9.56C-7 

9.77C-7 

9.99C-7 



1.290C-2 


1.02e-6 


2.621 

20.0 


3.029 

1.054C-2 

8,0e-7 

3.001 

1.063C-2 

8.09e-7 

2.946 

1.080e-2 

8.27C-7 

2.904 

1.098e-2 

8.45e-7 

2.865 

1.117e-2 

8.64e-7 

2.833 

1.137e-2 

8.82e-7 

2.805 

1.158e-2 

9,01e-7 

2.781 

1.178C-2 

9.19e-7 

2.760 

1.199C-2 

9.38e-7 

2.741 

1.220C-2 

9.56C-7 

2.725 

1.241C-2 

9.75e-7 

2.710 

1.262e-2 

9.92e-7 

2.697 

1.284e-2 

l.Ole-6 

2.684 

1.305C-2 

1.02e-6 

30.0 



3.261 

1.173C-2 

8.7e-7 

3.228 

1.178e-2 

8.8e-7 

3.171 

1.189e-2 

8.95e-7 

3.121 

1.201C-2 

9.06C-7 

3.072 

1.212e-2 

9.16C-7 

3.032 

1.225e-2 

9.29e-7 

2.998 

1.239C-2 

9.43C-7 

2.964 

1.253C-2 

9.57C-7 

2.937 

1.267e-2 

9.70C-7 

2.912 

1.281e-2 

9.84e-7 

2.887 

1.295C-2 

9,98e-7 

2.868 

1.310e-2 

l.Ole-6 

2.849 

1.324C-2 

1.02e-6 

2.83 

1.339e-2 

1.02e-6 


0.0 

0,04 

0.12 

0.2 

0.28 

0.36 

0.44 

0.52 

0.6 

0.68 

0.76 

0.84 

0.92 

1.0 


in-94 



pu zliq ztank dmdot 

20. 120. 120. .005 

************************************************************ 

zsm dsi zsi n nbar 

Jjc*********************************************************** 

132. .444 120. 45 4 

:tc ♦ 4( 4« 9<e 9k ♦ 4e 3|( 3<t ;(( 4e ♦ 4c 4( 9(< ♦ 9f( 9k ♦ 9|e 4( * 3f( 3k 9|( 9(t ** 4c )k ik ♦♦♦ 3ie ♦ ♦ 

ks acc rho vise roverd 

3k 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 9k 4c 4c 4c 4c 9k 9k 4c 4c 4c 9k 9k 9k 9k 4c 4c 9k 4c 9k 4c 9k 4c 4c 4c 4c 9k 4c 4c 4c 4c 4c 9k 9k 9k 9k 9k 4c 4c 4c 4c 9k 4c 4c 4c 4c 4c 4c 4c 4c 

.5 1. 4.339 0.817e-5 2. 

*«*****««:li*]|t%4c*4i****«***«ilt**]|t«#«*****3.*******)|iilti.*****il*4t*il<4t** 

tol nUm 

4 c 4 c 4 i 4 c 4 i 4 c 4 c 4c9k 4 c 4 c 9k 9k 9k 4c 4c 9k 4> 4c 9k 4c 4c 4c 4c 4c 9k 4c 4c4c 9k 4c 4c 9k 9k 4c 4c 4i 4c 4c 4c 9k 4c 9k 4c 4c 4c 4c 9k 4c 4c 4c 4c 4c 4c 4c 4c 9k 4c 4c 9k 

.001 50 

3k 4c 9k 4c 4c 9k 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 9k 4c 4c 4c 4c 4c 4c 4c 4( 4c 4c 4c 9k 4c 9k 4c 4c 4c 4c 9k 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4< 

opn opt omdin omds ovel oedas 

3k 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 9k 4c 4c 4c 4c 4c 9k 4c 4c 4c 9k 9k 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 9k 4c 4c 4c 4c 4( 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 

0 . 0 . 0 . 0 . 0 . 0 . 

20 PSIA ULLAGE, 0.2 LBM/SEC, 1 G. FULL TANK, 45 ORF 
DISTANCE FROM SPRAY TUBE INLET (IN) 

PRESSURE (PSIA) 

FLOW RATE THROUGH SPRAY TUBE (LBM/SEQ 
INJECTION FLOW RATE (LBM/SEQ 
INJECTION VELOCITY (FT/SEQ 
ORinCECDA(IN2) 

8 

4 .001901 
8 .002224 
12 .002534 
16 .002821 
20 .003079 
25 .003629 
30 .004122 
45 .004647 

3.5.2 lDig.gratsdZgrQ‘8 TVS Modd 

4t4c4(4c4c4c4c4c4c4c4c4c4c4c4c4c4c4c4c4c4c4c4c4c4i4>4c4c4c4j4c4c4c4(4c4c4c4c4c4c4c4c4c4c4c4c4c4c4c4c4c4c4c4c4c4c4c4c4c4c 

xd xchar he xcond prtsp outp input 

4c4c4c4c9k4c4c4c4c9k4c4c4c4c4c4c4c4c4c9k4c9k4c4c4(4c9k9k4t9k9k9k9k4c4c4c9k9k4c9k9k9k9k9k9k9k9k9k9k4c9k9k4c4c9k9k9k4c9k9k (\nfn 

1.0 0.25 0.0 1.0 13200.0 1.0 

4c 4c 4c 4t 4c 4c 4c 4c 4c 4c 4> 4c 4> 4c 4c 4c 3k 4> 4c 4c 4c 4c 4c 4c 4c 4c 4c 4( 4c 4( 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4( 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4< 4c 4c 4c 4c 4c 4c 4c 



pui tui pli twi twli full xll 

% 3ic 9|c 9(e )|c 9k 4< )|( >k 9k 3(c lie ak 9(( 9tc ♦ 4c 3k licik Ik * * * * ♦ 9k 4c 4t * 4c Jk 4c 3<C 9|( 

20.0 38.431 20.0 38.841 38.841 10.0 5.375 

4 e 4 c 4 e 4 c 4 e 4 e 4c 4c 4c 4c 4c 4c 4( 4c 4( 4( 4c 4( 4c 4( 4c 4c 4c 4c 4c 4c 4c 4c 4c 4( 4c 4c 4« 4( 4c 4c 4c 4c 4c 4c 4( 4c 4c 4c 4c 4« 4c 4c 4c 4c 4c 4c 4c 4c 4c 4( 4c 4( 4< 4( 4( 4c 4c 4c 4c 4c 4( 4c 4c 


dtank hcyl hbulk tkw dsb dl 


d2 


4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4« 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4( 4c 4c 4c 4c 4c 4( 4c 4c 4( 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 

120.0 60.0 30.0 0.5 0.25 1.902 1.402 

9k 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4t 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4( 4c 

mdvent mdsi dthex qflux g hwliq 

4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 

0.0052 0.001 2.2167 0.25 l.Oe-6 25.0 

pmin pmax delt2 ipmt2 iplot2 

4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 


20.0 21.0 1.0 15 


15 


4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 

fintim deltl xdeltl ipmtl iplotl nline 

4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 


12500.0 0.1 0.1 10 10 


40 


4c 9k 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4( 4c 4( 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 


opu 


otu 


omu opl otl oml 


4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4* 4c 4c 4t 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4t 4c 4c 4c 4c 4c 4c 4c 

1.0 0.0 0.0 1.0 0.0 0.0 

4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 9k 4c 4c 4c 4c 9k 9k 9k 4c 9k 9k 4c 9k 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 

otw omwl omdv omdsu omddu omdbw 

9k4c4c9k9k9k4c9k4c4c4c4c4c9k4c4c9k4c4c4c4c4c4c4c4c4c4i4c4c4c4c4c4c4c4c4c4c4c9k4c4c4c4c4c4c4c4c4c4c4i4c4c4c4c4c4c4c4c4c4c 

0.0 0.0 1.0 0.0 0.0 0.0 

4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 9k 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 9k 4c 4c 4c 4c 4c 9k 9k 4c 4c 4c 9k 4c 4c 4c 4c 4c 4c 4c 9k 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 

omdlu omdul omds onpump odppmp 

4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c4c 4c 4c 4c 4c 4c 4c4c4c 4c 4c 4c 4c 4c 4c4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 

0.0 0.0 1.0 1.0 1.0 


TIME (SEC) 

ULLAGE PRESSURE (PSIA) 

ULLAGE TEMPERATURE (R) 

ULLAGE MASS (LBM) 

BULK UQUID PRESSURE (PSIA) 

BULK LIQUro TEMPERATURE (R) 

BULK UQUID MASS (LBM) 

WALL TEMPERATURE (ULLAGE SIDE) (R) 
WALL UQUID MASS (LBM) 

OVERBOARD VENT FLOW RATE (LBM/SEQ 
ULLAGE SPRAY FLOW RATE (LBM/SEQ 
UQUID DROPLET BOILING RATE (LBM/SEQ 
WALL LIQUID BOILING RATE (LBM/SEQ 


r 


S' ; 


F 


m-96 



INTERFACIAL MASS-TRANSFER RATE (LBM/SEC) 
PUMP FLOW RATE (LBM/SEC) 

PUMP SPEED (RPM) 

PUMP PRESSURE RISE (PSI) 

ZERO-g TVS TRANSIENT PERFORMANCE 
(Og, 10% FULL, 0.25 BTU/HR-FT2, NO He) 


**j|c:|c********j|e*Jt«********************************************* pUITlp.f 

mdotd dpd npumpi npumpd xhp xn input 
*****************************************************’*‘****** ddtE 

0.3 0.500 0.0 3134.0 5.0 0.01 

******>►*♦**♦♦♦♦♦*♦****♦*♦***♦*******♦*********************** 

deltat effp 

************************************************************ 

0.2 0.65 

♦♦♦♦♦♦★★♦♦♦★♦♦♦♦♦♦♦★♦♦♦♦♦♦♦♦♦♦♦***************************** spray .f 

dsm zsm dsi zsi norf nbar input 

*♦******♦*****♦*♦♦*♦*♦*★********♦**★♦♦******♦*****♦♦*******♦ data 

1.18 132.0 0.444 120.0 45 4 

♦ ♦♦♦jIi******************************************************* 

ks cds roverd dmdot tol nlim 

♦*♦♦**♦***♦♦♦♦♦♦**♦**************************♦*************♦ 

1.56 0.8 2.0 0.005 0.001 200 

1 

45 0.0670 

24 p h cv cp rho beta k mu Saturated hydrogen 
properties 

24.845 1.021 -132.892 1.126 1.557 4.80827 .0058609 .04199 1.751e-5 

60.357 1.484 2.513 .00784 .0417855 .00719 .050e-5 

25. 1.073 -132.647 1.129 1.568 4.80377 .0059353 .04250 1.730e-5 

60.699 1.484 2.518 .00819 .0415883 .00721 .050e-5 

26. 1.462 -131.030 1.151 1.619 4.77434 .0061184 .04551 1.605e-5 

62.879 1.487 2.532 .01077 .0404077 .00739 .052e-5 

27. 1.950 -129.360 1.174 1.669 4.74429 .0063005 .04800 1.496c-5 

65.002 1.491 2.551 .01391 .0393821 .00756 .054e-5 

28. 2.553 -127.633 1.198 1.723 4.71359 .0064909 .05004 1.400e-5 

67.062 1.496 2.573 .01765 .0384993 .00775 .057e-5 

29. 3.288 -125.846 1.221 1.786 4.68219 .0067848 .05168 1.315e-5 

69.056 1.501 2.598 .02207 .0377489 .00793 .059e-5 

30. 4.170 -123.995 1.245 1.849 4.65003 .0070405 .05297 1.240e-5 

70.977 1.507 2.627 .02724 .0371221 .00813 .061e-5 

31. 5.217 -122.077 1.267 1.915 4.61705 .0073258 .05402 1.172e-5 

72.821 1.513 2.659 .03322 .0366118 .00834 .063e-5 


in-97 



32. 6.446 -120.090 1.289 1.985 4.58321 .0076421 .05489 1.112e-5 

74.584 1.520 2.695 .04008 .0362124 .00856 .066e-5 

33. 7.877 -118.029 1.310 2.051 4.54844 .0079185 .05555 1.056e-5 

76.261 1.527 2.734 .04790 .0359198 .00880 .068e-5 

34. 9.527 -115.893 1.329 2.125 4.51266 .0082741 .05601 1.006e-5 

77.848 1.535 2.778 .05675 .0357318 .00904 .070e-5 

35. 11.416 -113.678 1.348 2.195 4.47582 .0085840 .05631 .959e-5 

79.339 1.543 2.826 .06671 .0356477 .00929 .072e-5 

36. 13.561 -111.380 1.365 2.270 4.43782 .0089505 .05690 .917e-5 

80.729 1.551 2.879 .07787 .0356683 .00962 .075e-5 

37. 15.984 -108.997 1.380 2.360 4.39858 .0094522 .05748 .877e-5 

82.105 1.559 2.935 .09011 .0357690 .00998 .077e-5 

38. 18.694 -106.524 1.395 2.443 4.35801 .0099007 .05793 .840e-5 

83.256 1.567 2.998 .10395 .0360086 .01036 .079e-5 

39. 21.723 -103.956 1.408 2.532 4.31600 .0103939 .05824 .805e-5 

84.290 1.576 3.068 .11930 .0363705 .01076 .082e-5 

40. 25.089 -101.289 1.420 2.637 4.27243 .0110382 .05843 .773c-5 

85.199 1.584 3.146 .13629 .0368634 .01117 .084e-5 

41. 28.813 -98.517 1.431 2.743 4.22718 .0117045 .05851 .742e-5 

85.976 1.592 3.233 .15504 .0374994 .01159 .087e-5 

42. 32.915 -95.636 1.441 2.848 4.18010 .0123677 .05848 .713e-5 

86.614 1.601 3.331 .17568 .0382948 .01204 .089e-5 


43. 37.415 -92.637 1.451 2.969 4.13102 .0131689 .05835 .685e-5 

87.103 1.610 3.441 .19839 .0392710 .01251 .092e-5 


44 , 

42.334 

-89.513 1.459 

3.097 

4.07975 

.0140457 

.05811 

.659e-5 



87.431 1.619 

3.566 

.22335 

.0404560 

.01301 

.095e-5 

45. 

47.693 

-86.254 1.467 

3.242 

4.02608 

.0150803 

.05778 

.634C-5 



87.584 1.629 

3.709 

.25078 

.0418864 

.01354 

.097e-5 

50. 

81.838 

-67.493 1.507 

4.307 

3.71079 

.0235331 

.05459 

.519C-5 


85.043 1.693 4.919 .43621 .0551860 .01694 .112e-5 


55. 130.299 -42.122 1.564 7.528 3.25720 .0531428 .04960 .410e-5 

73.413 1.831 9.187 .76511 .1048799 .02422 .135e-5 


8 11 h 

cv 

cp rho beta k 

mu 

T Superheated hydrogen properties 

0. .01 .02 .03 

.04 .05 

.06 .07 .08 .09 .1 

600. 


1. 25.000 






60.765 

1.483 

2.511 

.00763 

.00721 

.050e-5 

25.000 

75.138 

1.481 

2.492 

.00617 

.00804 

.063C-5 

30.750 

89.438 

1.480 

2.483 

.00518 

.00915 

.074C-5 

36.500 

103.700 

1.480 

2.478 

.00447 

.01055 

.086e-5 

42.250 

117.938 

1.480 

2.474 

.00393 

.01202 

.097e-5 

48.000 

132.160 

1.480 

2.473 

.00351 

.01347 

.108e-5 

53.750 

146.372 

1.480 

2.471 

.00317 

.01477 

.117e-5 

59.500 

160.577 

1.480 

2.470 

.00297 

.01607 

.128e-5 

65.250 

174.782 

1.481 

2.470 

.00265 

.01737 

.137e-5 

71.000 

187.139 

1.484 

2.472 

.00248 

.01850 

.146e-5 

76.750 


ra-98 



203.223 

1.489 

2.477 

.00228 

.01996 

.156e-5 

82.500 

10. 30.806 






78.249 

1.537 

2.790 

.05926 

.00704 

,071e-5 

34.263 

93.566 

1.502 

2.649 

.04952 

.01030 

.082e-5 

39.920 

108.372 

1.493 

2.591 

.04274 

.01168 

.093e-5 

45.578 

122.930 

1.489 

2.558 

.03766 

.01306 

.103e-5 

51.235 

137.339 

1.487 

2.537 

.03370 

.01439 

.114C-5 

56.892 

151.651 

1.486 

2.523 

.03051 

.01564 

,124c-5 

62.550 

165.892 

1.485 

2.513 

.02788 

.01692 

.133C-5 

68.207 

180.088 

1.486 

2.506 

.02570 

.01819 

.142e-5 

73.865 

194.258 

1.489 

2.503 

.02382 

.01944 

.151e-5 

79.522 

208.423 

1.493 

2.504 

,02220 

.02071 

.160e-5 

85.179 

222.605 

1.502 

2.509 

.02080 

.02198 

.169e-5 

90.837 

15. 36.608 






81.613 

1.556 

2.912 

.08508 

.00984 

.076e-5 

36.608 

97.298 

1.511 

2.722 

.07143 

.01107 

.087C-5 

42.242 

112.381 

1.499 

2.640 

.06178 

.01239 

.098e-5 

47.876 

127.120 

1.493 

2.595 

.05459 

.01374 

.108e-5 

53.510 

141.657 

1.490 

2.566 

.04898 

.01500 

.118e-5 

59.144 

156.056 

1.489 

2.546 

.04445 

.01625 

.128e-5 

64.778 

170.360 

1.488 

2.532 

.04072 

.01751 

.137e-5 

70.412 

184.599 

1.489 

2.524 

.03759 

.01876 

.147e-5 

76.045 

198.803 

1.492 

2.519 

.03491 

.02001 

.155e-5 

81.679 

212.996 

1.498 

2.519 

.03260 

.02126 

.164e-5 

87.313 

227.201 

1.507 

2.524 

.03057 

.02252 

.172e-5 

92.947 

20 38.444 






83.731 

1.571 

3.029 

.11058 

.01054 

.080e-5 

38.444 

99.896 

1.519 

2.790 

.09269 

.01169 

.091C-5 

44.060 

115.240 

1.504 

2.689 

.08029 

.01297 

.102C-5 

49.675 

130.163 

1.498 

2.631 

.07106 

.01427 

,112e-5 

55.291 

144.826 

1.494 

2.594 

.06384 

.01550 

.1210-5 

60.906 

159.318 

1.491 

2.568 

.05801 

.01674 

.1310-5 

66.522 

173,683 

1.490 

2.552 

.05327 

.01799 

.1400-5 

72.137 

187.979 

1.492 

2.541 

.04921 

.01923 

.1490-5 

77.753 

202.224 

1.495 

2.534 

.04574 

.02047 

.1580-5 

83.368 

216.452 

1.502 

2.533 

.04273 

.02171 

.1670-5 

88.984 

230.687 

1.512 

2.538 

.04010 

.02296 

.1750-5 

94.600 

25 39.975 






85.177 

1.584 

3.144 

.13584 

.01116 

.0840-5 

39.975 

101.806 

1.526 

2.857 

.11368 

.01224 

.0940-5 

45.575 

117.424 

1.509 

2.734 

.09846 

.01348 

.1050-5 

51.176 

132.526 

1.501 

2.665 

.08719 

.01473 

.1150-5 

56.776 

147.319 

1.497 

2.620 

.07839 

,01594 

.1250-5 

62.376 

161.903 

1.493 

2.590 

.07131 

.01717 

.1340-5 

67.976 

176.339 

1.493 

2.570 

.06554 

.01841 

.1430-5 

73.577 


in-99 


190.690 

1.494 

2.556 

.06058 

.01964 

.152e-5 

79.177 

204.980 

1.498 

2.548 

.05634 

.02087 

.161e-5 

84.177 

219.242 

1.505 

2.546 

.05269 

.02210 

.169e-5 

90.377 

233.514 

1.517 

2.550 

.04952 

.02335 

.177e-5 

95.978 

30. 41.299 






86.182 

1.595 

3.261 

.16101 

.01173 

.087e-5 

41,299 

103.265 

1.533 

2.924 

.13441 

.01274 

.098e-5 

46.886 

119.156 

1.514 

2.779 

.11635 

.01393 

.108e-5 

52.473 

134.439 

1.505 

2.697 

.10303 

.01523 

.117e-5 

58.060 

149.354 

1.499 

2.646 

.09273 

.01634 

.127e-5 

63.647 

164.030 

1.495 

2.610 

.08442 

.01756 

.136C-5 

69,234 

178.543 

1.495 

2.587 

.07756 

.01878 

.146e-5 

74,821 

192,947 

1.496 

2.571 

.07178 

.02000 

.155e-5 

80.408 

207.282 

1.500 

2.562 

.06685 

.02122 

,163e-5 

85.995 

221.586 

1.509 

2.559 

.06257 

.02245 

.171e-5 

91.582 

235.893 

1.521 

2.563 

.05881 

.02369 

.179e-5 

97.169 

35. 42.475 






86.865 

1.605 

3.381 

.18169 

.01226 

.090e-5 

42.475 

104.432 

1.539 

2.987 

.15482 

.01319 

.lOOe-5 

48.050 

120.569 

1.518 

2.822 

.13404 

.01434 

.llOe-5 

53.626 

136.023 

1.508 

2.730 

.11876 

.01551 

.120e-5 

59.201 

151.065 

1.501 

2.670 

.10690 

.01671 

.129e-5 

64.776 

165.833 

1.498 

2.630 

.09737 

.01791 

.139e-5 

70.351 

180.418 

1.497 

2.604 

.08953 

.01912 

.147e-5 

75.927 

194.879 

1,499 

2.586 

.08291 

.02033 

.156e-5 

81.502 

209.264 

1.503 

2.576 

.07724 

.02155 

.165e-5 

87.077 

223.610 

1.512 

2.572 

.07231 

.02277 

.173e-5 

92.652 

237.956 

1.524 

2.575 

.06799 

.02399 

.181e-5 

98.228 

40. 43.536 






87.300 

1.615 

3.506 

.21148 

.01278 

.093e-5 

43.536 

105.311 

1.545 

3.056 

.17544 

.01362 

.103e-5 

49.101 

121.728 

1.523 

2.866 

.15163 

.01472 

.113e-5 

54.665 

137.369 

1.512 

2.761 

.13427 

.01587 

.122e-5 

60.230 

152.534 

1.503 

2.694 

.12088 

.01705 

,132c-5 

65.795 

167,379 

1.500 

2.651 

.11027 

.01823 

.140e-5 

71.359 

182.041 

1.499 

2.621 

.10139 

.01943 

.149e-5 

76.924 

196.562 

1.500 

2.601 

.09390 

.02063 

.158e-5 

82.488 

210.999 

1.506 

2.589 

.08750 

.02185 

.166C-5 

88.053 

225.389 

1.515 

2.584 

.08193 

.02306 

.175C-5 

93.618 

239.769 

1.528 

2.586 

.07706 

.02428 

.183e-5 

99.182 



